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The Hubbard model with finite on-site repulsion U is studied via the functional-integral for- 
10 ' mulation of the four-slave-boson approach by Kotliar and Ruckenstein. It is shown that a correct 

0^ ' treatment of the continuum imaginary time limit (which is required by the very definition of the 

functional integral) modifies the free energy when fluctuation (1/A'') corrections beyond mean-field 
• are considered, thus removing the inconsistencies originating from the incorrect handling of this 

I pathologic limit so far performed in the literature. In particular, our treatment correctly restores 

D . the decrease of the average number of doubly occupied sites for increasing U. Our analysis requires 

pL^ ' us to suitably interpret the Kotliar and Ruckenstein choice for the bosonic hopping operator and to 

, abandon the commonly used normal-ordering prescription, in order to obtain meaningful fiuctuation 

■ corrections. In this way we recover the exact solution at (7 = not only at the mean-field level but 

also at the next order in In addition, we consider alternative choices for the bosonic hopping 

^SJ , operator and test them numerically for a simple two-site model for which the exact solution is read- 

^ ' ily available for any U. We also discuss how the 1/N expansion can be formally generalized to the 

OA ' four-slave-boson approach, and provide a simplified prescription to obtain the additional terms in 

. the free energy which result at the order 1/7V from the correct continuum limit. 

^ ■ PACS numbers : 71.30.-fh, 05.30.-d, 71.10.+X 

o ■ 
g: 

I i The Hubbard model (with its variations) is often adopted to represent the essential correlations among electrons 
■ O " (holes) in a lattice, especially in the intermediate- and strong-coupling regimes (relevant to heavy-fermion and high-Tc 
C materials). In these regimes, the use of conventional many-body techniques (that rely on expansions in terms of the 
^ on-site repulsion U) becomes questionable and alternative non-perturbative methods are required.^ Specifically, in the 
strong-coupling limit the need for projecting out the electronic configurations with double occupancy at any given site 
acquires the role of a "constraint" which has to be enforced for a correct description of correlations. Special methods 
have accordingly been developed to satisfy the local constraints in the strong-coupling limit .□ 
5^1 ■ In particular, slave-boson representations have been introduced to allow for field-theoretical treatments of strong 
' correlations with constraints, amounting to a projection method that employs auxiliary (or "slave") bosonic particles.u 
While in the original version slave bosons were referring only to otnpty and doubly occupied states at any given lattice 
site, the method was later extended by Kotliar and RuckensteinD (KR) who assigned slave bosons also to the singly- 
occupied states. The KR four-slave-boson method maps the physical fermion (destruction) operator fi^^ with spin a 
at site i onto the product of a (pseudo) fermion fi^^ with the same spin and of a bosonic operator Zi^^ (see below), 
and is especially suited to deal with the Hubbard model for any U. The KR method is also particularly appealing for 
the treatment of magnetic problems, which can be approached already at the mean-field-level owing to the presence 
of the single-occupancy bosons. For these reasons the KR method has beeru adopted to treat several problems, 
both at its mean-field levelEl and with the inclusion of fluctuation corrections Jj by relying on a functional-integral 
formulation which is ideally suited to implement the local constraints for slave bosons .0 Specifically, it has been found 
that the KR mean-field solution gives remarkable agreement with more elaborate Monte Carlo results over a wide 
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range of interactions and particle densities .0 This success is not entirely surprising since the KR (paramagnetic) mearu 
field solution has been modeled after the Gutzwiller results and thus is identical to tbe Gutzwiller approximationa 
(the latter having been further shown recently to become exact in infinite dimensionsll3) . The KR finding that the 
Gutzwiller results can be recovered without explicit wave functions has then prompted the hope to improve the theory 
systematically by the inclusion of fluctuation corrections, using the KR saddle-point solution as a good starting point. 

The use of the KR four-slave-boson method was for us originally motivated by the interest in describing the magnetic 
properties of Cu02 layers with an itinerant approach in the limit of large repulsion Ud at Cu sites. Antiferromagnetic 
calculations at the mean-field level within a three-band Hubbard model with C/^ « oo (where the double occupancy 
boson. does.not enter) were actually quite promising.til However, our attempts to include fluctuations along the lines of 
RefsJj andQ met with inconsistencies, leading to an instability of our previous mean-field results. Further consideration 
of the single-band Hubbard model with finite U did not remove the unpleasant features of fluctuation calculations (for 
instance, the average number of doubly occupied sites away from "half filling" was found to unphysically increase for 
increasing U with the inclusion of fluctuations). We were thus unavoidably led to conclude that there was something 
systematically wrong in the standard procedure adopted in the literature to include fluctuation corrections within the 
KR method. |— ■ 

About at the same time Jolicoeur and Le Guilloiila signaled the occurrence of additional inconsistencies when 
including fluctuation corrections within the KR method at U — 0, and proposed to modify the bosonic hopping 
operator z mentioned above at each order in the loop expansion (without suggesting, however, any explicit form of z 
even at the one- loop order). 

On the other hand, it was soon clear that the inconsistencies found when including fluctuations were not limited to 
the KR four-slave-boson method (and thus to the presence of the operator z whose choice is to a large extent arbitrary) . 
In fact, Emery and coworkeraiS have shown that inconsistencies arise when including fluctuations even in the simpler 
U = oo one-slave-boson problem (where only the empty boson appears), which is free from the arbitrariness related 
to the operator z. This result then suggested that the origin of the problems was in the slave-boson approach itself, 
irrespective of its particular version. 

The solution to these problems was obtained eventually by deeper examination of the functional-integral formula- 
tion on which slave-boson methods rely. Speciflcally, it turned out that the continuum imaginary time limit of the 
functional-integral representation of the partition function had been incorrectly implemented beyond the mean-field 
level in the slave-boson literature, by performing the continuum limit in the action at the outset. By doing so, the 
bosonic commutators were not properly represented in the functional integral, a shortcoming which turns out to be 
crucial in the presence of (slave-) boson condensates. The peculiarities of the continuum imaginary time limit within 
a coherent-state functional-integral approach in the presence of a slave-rbpson condensate have been originally demon- 
strated in the context of the simpler U — oo one-slave-boson problem.Ej In the present paper we concentrate on the 
implications of taking the correct continuum limit with the four-slave-boson method, which prove to be nontrivial 
due tOptte presence of the bosonic operator z. A short preliminary version of the results presented here can be found 
in Ref.lla 

Taking the continuum time limit at the end of the calculation is actually a well-established procedure for Feynman's 
path integrals.tj On the other hand, for the functional integrals based on the coherent-state representation it has been 
common practice to consider the continuum limit at the outset, although warnings have been given that this procedure 
might lead to inconsistencies .liZI Taking the continuum limit at the outset considerably simplifies the calculations, by 
making the expression of the action more manageable and allowing for the use of standard Matsubara techniques 
developed for the diagrammatic Green's functions approach.li3 Conversely, keeping a finite imaginary time mesh 
until the end of the calculation requires one to reconsider the Matsubara techniques for a finite set of (fermionic. 
or/and bosonic) frequencies. Although this task has been avoided in the previous slave-boson literature, our workllil^Ej 
demonstrates unambiguously that, in the presence of bosonic condensates, it is necessary to preserve the discretized 
form of the action until the end of the calculation when resorting to a coherent-state functional-integral representation 
of the partition function. In fact, proper account of the discretized nature of the functional integral yields additional 
terms for the free energy when including fluctuations, which would otherwise be missed by taking the continuum limit 
at the outset. 

Concerning further the KR four-slave-boson method, the presence of several bosons and of the bosonic operator z 
makes unavoidablK-jthe derivation of the above additional terms more involved than for the one-slave-boson method 
considered in Refo. These additional terms turn out to be essential for both methods to heal the inconsistencies 
found in the literature when taking the continuum limit at the outset. In the present paper we derive in detail 
these additional terms for the four-slave-boson method with a generic form of z. In particular, we will show that 
the form of the bosonic operator z proposed by Kotliar and Ruckenstein (which, for the sake of a direct mapping 
onto a functional-integral formulation, has been invariably interpreted in the literature within a normal-ordering 
prescription and which will be referred to as zk r in the following) leads to incpasistencies when including fluctuations 
even by taking carefully the continuum limit at the end of the calculationJlS Nonetheless, we will also show that 
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these inconsistencies can be overcome by suitably relaxing the normal-oijdering prescription on the KR choice for the 
bosonic operator z. In this way, the preliminary results presented in Ref.t£l are extended and improved. In particular, 
this form of z without the normal-ordering prescription (which we shall refer to as zsq) turns out to reproduce the 
correct independent-particle solution at [/ = not only at the mean-field level but also with the inclusion of the 
corrections. (Numerical results will be considered in the zero-temperature limit only throughout this paper). In 
addition, we shall consider for comparison two alternative forms of z (different from zkr and zsq) which are also free 
from the inconsistencies occurring for zkr- Specifically, we will consider a "linearized" form zlin which preserves, 
by construction, the KR mean-field solution at "half filling" of the paramagnetic band for any value of U and also 
strongly suppresses the contribution of fluctuations to the free energy at J7 = (although not completely and at "half 
filling" only). Our finding that the results obtained for the ground-state energy alternatively with zgq and zl/at do 
not differ appreciably for large values of U (say, for U > 4t) even away from "half filling" hints further that the results 
obtained with the four-slave-boson method might not depend crucially on the choice of the operator z for practical 
purposes. 

Numerical calculations will be mainly carried out for a one-level two-site model, which avoids the complications 
due to the spatial structure and keeps those due to the imaginary time discretization of the functional integral which 
we are mostly interested in. Restriction to a two-site model will enable us to compare our numerical results, obtained 
within alternative approximations to the functional integral, with the exact solution which is readily available for any 
value of U and band filling. 

The plan of the paper is the following. Section |l| deals with the formal problems related to the correct functional- 
integral formulation of the four-slave-boson method over a discretized imaginary time mesh, and with the peculiar 
problems related to the choice of the bosonic factor z. The novel terms for the free energy, resulting at the Gaussian 
level from taking the continuum limit of the functional integral only at the end of the calculation, are derived 
in detail for any given form of the bosonic factor z and interpreted by the need of restoring the correct bosonic 
commutators. Section III presents our numerical calculations with the four-slave-boson method for the ground-state 
energy of a simple two-site model and compares them with the available exact results for any IJ and band filling. 



The results obtained with alternative forms of the bosonic operator z are discussed in this context. Section [V gives 
our conclusions. Finally, the Appendices contain additional mathematical arguments which are needed for making 
the material presented in the text self-contained. In particular. Appendix |A| discusses the XjN expansion for the 
four-slave-boson method, and Appendix |c| recovers within the "Cartesian" gauge the additional terms for the free 
energy (which were derived in the text within the "radial" gauge). The problem of obtaining the additional terms 
(due to a correct handling of the continuum limit of the functional integral ) also for the correlation functions (and 
not only for the free energy) is briefly discussed in this context. 



II. FUNCTIONAL-INTEGRAL FORMULATION OF THE FOUR-SLAVE-BOSON METHOD 

In this Section, we discuss the procedure to formulate the four-slave-boson method via the coherent-state functional- 
integral representation of the partition function, by keeping the discretized imaginary time mesh (which is required 
by the very definition of the functional integral) until the end of the calculation. Specifically, we will demonstrate at 
the Gaussian level that this correct procedure yields additional terms for the free energy, which were missed in the 
previous literature when the continuum limit of the functional integral was taken incorrectly in tbe. effective action 
at the outset. By doing so, we extend to the four-slave-boson method the results obtained in RefJlj for the simpler 
U = oo one-slave-boson problem. Mastering the formal apparatus with a discretized time mesh, however, is now 
unavoidably more involved. 

For the sake of definiteness, we consider the single-band Hubbard Hamiltonian 

i,A,fj i 

with on-site repulsion U, where cr is a spin label and A runs over the "star" of nearest neighbors to site i in a two- 
dimensional square lattice of M sites. Extension of the following arguments to more complex lattices with multi-band 
structures and to off-site repulsive terms should, in principle, be straightforward. 

As mentioned in the Introduction, an essential requirement on any method for an approximate solution of in the 
strong-coupling regime {U ^ t) is that double occupancy at a given site is properly treated. To this end, suppression 
of double occupancy in the limit U/t — > oo acquires the role of a "constraint" for the many-body problem, whose 
enforcement poses no problem for systems of small size but is difficult to implement for an infinite system. In this 
case, the introduction of auxiliary variables in the form of slave bosons, which keep track of the occupancy at any 
given lattice site, looks particularly convenient for enforcing the "constraint" . 
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A. THE FOUR-SLAVE-BOSON METHOD 



For any given U, Kotliar and RuckensteinS have mapped the physical fermion operator /j.o- in Eq. ([^) onto the 
product fi^a-Zi^a of a (pseudo) fermion fi^a- and of a bosonic operator Zi_^ (at any given site). In its simplest version 
Zi rr has the form 



- J ^.^J 



+ elsi,<^ ' (2) 

where the bosonic destruction operators e^, Si^a, and di refer to empty, singly, and doubly occupied states (in the 
order) at site i. To establish a one-to-one correspondence between the original Fock space and the enlarged one which 
contains also the bosonic states, the following constraints must be satisfied: 

4^^ + 4,aSt.<y + ele, = 1 (3a) 

fUk^ - 4,.*^,- ~ 4d^ = 0. (3b) 

However, enforcement of the constraints (|^) allows for alternative choices of Zi^^ different from (^. This observation 
has been exploited by Kotliar and Ruckenstein who have accordingly modified the form (^ as follows 

Z'i^a — —(jRi.fjdi "t- G,^ u Si^(j , (4) 

with the requirement that the subsidiary operator Ri^a- acts as the unit operator in the subspace specified by the 
constraints (H). In particular, Kotliar and Ruckenstein have considered the nonlinear form 

i?...= =: : (5a) 



where : O : denotes the normal ordering of the operator O. This expedient has enabled Kotliar and Ruckenstein 
to recover the Gutzwiller solution for a paramagnetic band in a straightforward way, by considering the mean-field 
solution whereby each bosonic operator is replaced by its average value. Other forms of Ri^a different from the 
KR choice (^^ are, however, possible in principle. We shall exploit this freedom in the following, and explore also 
alternative expressions containing the bosonic number operators e|ej, sj^s^.o-, and djdi. 



B. FUNCTIONAL-INTEGRAL FORMULATION 



Functional integrals provide an ideal framework to enforce constraints like (|^). This is achieved by introducing a 
Lagrange multiplier for each constraiat (say, A" and A^^ for the constraints (^), in the order). The functional-integral 
formulation based on coherent statesE-O rests on breaking up the (imagiaary time) interval (0,/3) into M steps (where 
(3 = \/{kBT) is the inverse temperature and Af ^ oo eventually) ,E3 and yields the following expression for the 
grand-canonical partition function: 



M-l 



Z - Jrn^ f HdX'} [lldXl, j Yl d'e^.n'fd,,^, 

•' i \ a ) m=0 

X ^JJd^Si,o.,TOd/j,o-,md/i,<T,m^ exp{-S'Af}, (6) 

where m labels the imaginary time steps, (e, Sct,^) are complex boson fields, / and / are Grassmann variables, and 
Sm is the discretized action given by 

M~\ 

Sm^SJ:J: {w}'J + + - A^) . (7) 

m— i 

Here, 5 — P/M is the elementary time step. 
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i.m 



(8) 



(/I being the chemical potential) anc 



W^!™ - J<™ - (1 - ^5 A,") e,,™-i : 



1 



Si.a,m ~ ( 1 ^ '5 A^ + <5 _ J Si „ rn_i 



are functions of the Grassmann variables and of the boson fields, respectively, while 



(FB) 



(9) 



(10) 



is the mixed fermion-boson term which derives from the kinetic term in the Hamiltonian (|l|) and thus depends explicitly 
on the form of the bosonic operator Zi^^. [In the above expressions, the label a takes the two values ±1 corresponding 
to spin i. A generalized label S which can take 2N values (with A'' positive integer) will be introduced in Appendix 
^ as a formal tool to generate a l/N expansion for the four-slave-boson method]. 

In the expression (p^, Zi^„{m,m ~ 1) stands for the matrix element of the operator Zi^a- between coherent states 
labeled by imaginary times = rnd (on the left) and Tm-i — Tm — 5 (on the right). When Zi.o- is taken of the general 
form (U), this matrix element becomes 



Zj_o-(m, m') ^ [e 



Ri^a{m, m!) 



ill) 



where Ri^cj{'m, m') is the matrix element of the subsidiary operator Ri u between coherent states at times m and m' . 
In particular, when the normal-ordered form (|5a| ) is considered, the above matrix element Ri^n{m,'m') acquires the 
simple form: 



i?5^(m,m')- 



1 



1 



(12) 



This form with m = m' has been constantly used in the previous literature.0 In the more general case when the 
operator Ri^^ is not explicitly in normal-ordered form, it should be expressed as a sum of normal-ordered terms by 
suitably commuting the bosonic operators. In practice, for an operator like ( pb| ) this procedure can be implemented 
in the context of a 1/A^ expansion, as it will be shown in subsection II E . 

As remarked in the Introduction, the functional-integral representation of the paxtition function requires one to 
keep, in principle, the discretized imaginary time mesh until the end of the calculation.EZi However, it has been common 
practice in the slave-boson literature to violate this requirement by taking the continuum [5 — *■ 0) limit of the action 
(0) " (0) the outset, thus effectively transferring the Af — > oo limit in Eq. (^) under the integral sign. This 
procedure removes the distinction between the two time labels m and m! in Eqs. (|l^) and (p^). These seemingly 
formal considerations acquire relevance from the occurrence of the unphysical results obtained in the continuum 
limit already at the Gaussian level, as shown in the next Section. 

Specifically, we argue that taking the continuum limit at the outset is pathologic for the slavtirboson problems 
owing to the presence of the slave-boson condensate, as it was already discussed in detail in Refo for the U — oo 
one-slave-boson problem. Accordingly, the M — > oo limit has to be properly taken in Eq. (^ only at the end of 
the calculation. Otherwise, wrong results for the free energy (and derived quantities) are obtained when including 
fluctuation corrections beyond the saddle-point of the functional integral. As it will be clear from the following 
analysis, these problems stem from a failure of the functional integral to account for the bosonic commutation rules 
when the continuum limit is taken at the outset. 



C. GAUSSIAN FLUCTUATIONS FOR THE FREE ENERGY WITH A DISCRETIZED TIME MESH 

The complex boson fields at each discretized imaginary time Tm entering the functional integral (^ can be rep- 
resented, alternatively, either via their amplitude and phase or via their real and imaginary parts. The two repre- 
sentations correspond to the so-called "radial" and "Cartesian" gauges, respectively. In this Section, we adopt the 
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"radial" gauge to evaluate on the same footing both the conventional "continuum" contribution (which was regularly 
considered in the previous literature) and the novel contribution resulting from a careful handling of the continuum 
limit of the functional integral (which we shall name "contribution from infinity " from the range of Matsubara 
frequencies it originates from) . Accordingly, each stage of the calculation reported below will be free from infrared 
singularities which plague instead the "Cartesian" treatment. However, we shall show in Appendix ^ that use of the 
"Cartesian" gauge makes the specific calculation of the "contribution from infinity " somewhat simpler than in the 
"radial" gauge (at least at the Gaussian level we consider) . Comparison of the results obtained for the "contribution 
from infinity " in the two gauges can also serve as a test of the validity of the 1 /N expansion for the four-slave-boson 
method presented in Appendix ^ 

The boson fields in Eq. (0) are represented via their amplitude and phase, by setting 




e^,„ exp{i^|^2^} 



(13) 



dj,„exp{i(^m. 



It is then convenient to transform the Grassmann fields therein as follows: 



(14) 



Entering Eqs. ( p^ ) and (14) in the action (|^) - ( pO| ) makes the boson phases to appear only in certain combinations. 
It is then convenient to introduce the following alternative variables 



A° — A" — 



As.) 



-y v'^""^) 



(15) 



In this way, the three phases (fi'f^ and ^^^^ (with a 



±1) are eliminated in favor of the new variables Af^ and 



a nn which Can be considered as time-dependent Lagrange multipliers. The need of keeping in the functional 
integral an additional independent variable, like gi^m of Eq. (p^, has been already pointed out by Jolicoeur and Le 
Guillou (Ref.Ei) [although in the continuum case], by arguing tliRt (p^ ^ Ccinnot be chosGii to make Qi^rn 

vanishing for 

each m. For later convenience, we introduce also the variables 



\a _ \b 

i.rn—1 i,a.m—l 



(e) 
T i^rn — 

- ip. ' 



l) 



(16) 



which, however, are not linearly independent because 

A/-1 M-1 



0. 



m— m— 

To simplify the notation somewhat, we further introduce the dictionary 

= a^a(5)^;,(5) 



e <-!■ 



Si ^ a 



(2) 
(3) 



a(3), 



(17) 



(18) 



with the understanding that the appropriate site and time (or wave vector and frequency) indices (or arguments) will 
be indicated whenever necessary. In addition, the notation a will stand for the set { a(");a = 1, . . . , 8}, & for the set 
{6('^);/3= 1,...,5}, and A for the set {A^M^ 1,...,3}. 

With the above notation, the terms (|1)-(M) of the action (Q) can be cast in the form: 



{F) 



fi,cr,m Q (A^.o- m — 1 , Aj _ fJ-) fi,tT,m—l 



(19) 
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0=1 



(20) 



and 



(FB) 



(21) 



Here, q"^ , , and p'^ (p"^) arc functions of the bosonic and A variables, defined respectively by 



-''(\>,m-i > >H,a - m) = [l - S{Xl^ - /i)] exp {-(5A-_^_^_i } , 



(22) 
(23) 



h " (Ai,„^-l , Ai, g 



.1 ) = [1 - <5(A," - exp {-5 (A?^„_i ~ )} , 



(24) 



and 



/i'(A.,™-l , A., 5..™ - 9^,m-l ) = [1 - <5(A^ - ^ A^, J - JC/] 



exp |-5 (A^^„„i - ^ A^_„^„„i ) - (5,, 



m 5i , m — 1 ) r ; 



^"(0*,™, a^,™-! ) = exp{-5A^„_i I^K^^m, "I'm , "!'m ' 



(25) 



(26) 



with cr = — a and 



1,171 

Jd) _ 



-1 exp{-(5Af,„_i} 



.1 exp{-(5 (A^^„_i - A,"^ )} 



n; ™ = d^^„A,ra-i exp{-(5 (Af^„_i - A° ) - (5^,™ - )} 



(27) 



[ p'^ is obtained from p"' by interchanging 6,-'^,^ ^ ^I'^n-i (/^ = l^'^jS) and then letting g — > —17 ]. In Eg. ([2^) 
it is understood that the subsidiary Junction TZ is formally obtained from the subsidiary operator Ri o- of Eq.(^) by 
replacing each bosonic occupation number therein with the corresponding c- number argument ( |2^ ) . This prescription 
holds provided the subsidiary operator R is written as the normal ordering of a function of the bosonic occupation 
numbers. 

Equations (pj|)-(27) are still exact. We now proceed and expand the action in powers of the fluctuating flelds about 
a chosen mean-field solution, by setting for the bosonic and static A variables 



i.m 

Af ) = A[/) + A 



(0 



(/3=l,---5), 
(/ = !,... 3). 



(28) 



The choice (^) entails a homogeneous (i.e., site independent) mean-field solution, which will restrict us in the following 
to the paramagnetic (or else, to the ferromagnetic) case. It also implies that the mean- field solution is static (i.e., tim|e-| 
independent ).LJ Accordingly, the dynamic variables ( p^ ) have, by their definition, vanishing mean-field contribution.E^ 
We next introduce the space and time Fourier transforms of the fluctuating fields, by setting 



BZ M-l 
q 1^=0 



(29) 
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BZ M-1 



\(0 ^ \(0 

i.m i.m 



and 



q p—1 



BZ 



^(0^^g.q.R.^(0(g^O) (/ = !,• ••,3). 



(30) 



(31) 



Here, Ri is the lattice vector associated to site i, q is a wave vector restricted to the first Brihouin zone (BZ), r,„ = m5 
(to = 0, • • • , M — 1), = 2-Kvl (i is a bosonic Matsubara frequency, and A^'^(q, 0) = \^^\q, uj^ = 0) by our definition. 
A similar transformation holds for the Grassmann fields: 



BZ M-l 



(32) 



k s=0 



where now uj^ — 27r(s + l/2)//3 is a fermionic Matsubara fpnuency. Note that the range of the Matsubara frequencies 
remains hounded while keeping the discretized time mesh.Ej 

Expanding the action up to quadratic order in the fluctuating fields requires one to consider the first and second 
derivatives of the functions (|2^)-(|2^) with respect to their arguments. By our convention, each argument may represent 
either a single variable or a set of variables (in the latter case we shall underline the corresponding symbol) . To keep 
the notation compact, we then introduce the following short-hand notation for the required derivatives. We set: 



9o;0 = (f{x,y) 



1l;0 = 

(T 

%;1 = 

9(l,l);0 

90;(1,1) 



dx 



x=0,y=\''^„-fj. 



dy 



x=0,y=\lg-ti 



d^<f{x,y) 



d'^<f{x,y) 



dy^ 
d'^q"{x,y) 



x=0.,y=\lg-i-i 



dxdy 



(33) 



p dh^jX, Y,Z) 
. _dh^iX,Y,Z) 

"■0;0;1 ^ 



X=0, Y=X„, Z=0 



X=0, Y=Xg , Z=0 



dZ 

,0 _d^h^{X,Y,Z) 



,/3 „ h^{X, Y,Z) 

^/;0;l ~ 



X=0,Y=Xr,,Z=0 



dXO)dZ 



X=O.Y=X„,Z=0 



,/3 _ d^h^jX, Y,Z) 



92 Y, Z) 



^0;0;(1,1) 



X=0,Y=Xr.,Z=0 



x=o, ii=Aq, z=o 



(34) 



and so on; 
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P0;0 - P^i^^ — )U=(^o'2).r=(6o'0) 



Pa-Q 



Pa:ct' 



X = {b„.0).Y = ib„,Q) 



d^P^iX, Y) 



_ d^p^{X,Y) 



X = iba,0),Y={bo,0) 



2L={b„,Q),Y=ib„.0) 



(35) 



and so on. Similar conventions hold for p'^. In particular, the suffices X = (607O) ^-nd Y_ = (feg, 0) in Eqs.(^5|) signify 
that each bosonic field bf^ (/3 = 1,- • - ,5) is replaced by the associated mean- field value while each dynamic 

Mm ^^I'i = 1; • ■ • ! 3) is replaced by zero. 
Expanding at this point the action (0) up to quadratic order in the fluctuating bosonic fields, we obtain the following 



contributions: 



Here: 



= Kij: br - 1) - E ^^Ub'i'' + ^'') ^ub'i^' 

13=1 a- 

is the constant (mean-field) value of the bosonic part (||) of the action; 

BZ M-l 

4°^ =EE E y<rik,0Js) e"^=^G^(fc,c.,)-V.(fc,c..) 

where 



k a s=0 



(e- 



-iuJs 5 



1) 



is the single-particle mean-field fermionic Green's function corresponding to the band eigenvalue 

ea{k) = A^_o-M + ^Po;oPo;o7(fe) 

with 

7(fc) = E*5xp{-jfc • A} ; 



(36) 



(37) 



(38) 



(39) 



(40) 



(41) 



where 



BZ M-l 8 BZ M-l 



4^^ -E E EEE E 7^(fc'-^) e-^£:«(g, ..;/.,.!«) 



q I/— a— 1 h (7 s—0 



X a'"^ (q, cji.) /<^(fc - q, uj^ - t^i.) 



fW(q,c^.;fc,aH = (l-5,^oE<'+5) 



(=1 



+ 7(fc-g) P^;o(Po;a+ '^^pS^o) 



(42) 



(43) 
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5^^j (= 1 when i = j and otherwise ) being the Kronecker delta function; 

BZ M-1 8 

^^'^ = E E E "^"^ («' '^-) ^-1"' -'^-) 

q 1^=0 a,a' = l 



where 



B(g,a;.|a,a') = C' E<4 (1 ' ^~'^"'Kfi-fi) 



0=1 



1 (1 + e-"^-^) 



/3=1 



(' = 1 



('=1 



/3' = 1 



1=1 

1 ^ (-3 3 

~ ^ E ^0 ^ E ^a,5+l E "^^,5+;' 
[3=1 I (=1 i' = l 

^ ^1^,0 )'''f;,i');0;0 + ^1^,0 ''■0;(;,i');0 

+ E<5+ie,5(e^'^'''-l)/iCo;l 

+'5^5E^5,5+;'(e-^'^'''-l)/i(^;0;l 

;'=i 

+^5^,5(6''^*'' - l)(e--''^ - l)<o;(l,l) 



(44) 



(45) 



where 



BZ M-1 



4'^ = EE E «^"^(^.-^) 

q 1^=0 a,a' = l 



BZ M-1 

EEE7.(fe. 

_ k (7 s— 



w., e 



x£^'^\q, ujv] k, a\a, a')fa{k, uis) 



f (2)(q, u;.- k, a\a, a') = - <5^„,((5^,6 -5^+1 + -^a.T S^,-i) 
x^((l- ^^o)9ri,i);0+^^^o<(i,i)) 

3 3 



t 

^2 



^' = 1 



7(«) P0;0(P(a,a');0+ P0;(a,a')+ + Pa';a e ") 



(46) 
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(47) 



Note that in the action (|3^) we have not considered the linear term in the bosonic field (which originates from the 
purely bosonic contribution (||)), since this term cancels at self-consistency by definition. lj In addition, only pairs 
of fermionic (Grassmann) variables with the same k and tOg have been retained in the contribution (|4^), because all 
other pairs would give a vanishing contribution to the following Gaussian calculation. 

The fermionic variables are eliminated at this point by performing the functional integration over the action Sp"^ + 

(2) 

Sp which is quadratic in the fermionic variables. Exponentiating the resulting expression and expanding the 



exponent up to second order in the fluctuating bosonic fields, we eventually obtain the effective action 
where S'(o) and S'^^^ are given by ( ^T\ ) and (^), respectively, and 

BZ M-l 
BZ M-1 8 

^^'^=EE E a^"\q,u;.)C{q,u;,\a,a')a^'^'H-q,-oj,), 

q L/—0 a.a' — l 



(48) 

(49) 
(50) 



with 



^ BZ 

C{q, uj^\a, a')= ^ ^|/M(e^(fc))f (q, k, cr|a, a') 

k a 

-]-£'^^\q,uj^\k,a\a)I{M {£a{k),e„{k - q);v) 



<£(i)(-g,-c^,;fc-9,fT|a')} 



(51) 

[cf. Eqs. (^9[), (|43|), and ([17[)]. In expression (|5l|) we have introduced the Fermi function /M(e) and the fermionic 
polarization function TIm{£, £']i^)i which generalize to the case of the discretized time mesh considered in this paper 
the familiar functions of the Matsubara Green's functions theory, the latter being recovered in the continuum [6 0) 
limit. The expressions of /m(£) and nM(e, i^) are reported in Appendix 

The free energy is obtained by performing the Gaussian integration over the bosonic variables { a*^"^ (q, w,y) ; a — 
1, • • • ,8} with action Seff given by Eq. (^8|). One obtains (per lattice site): 



F — _Fo + _Fi 



with 



and 



Fq = + 5' 



(0) 



(52) 
(53) 



BZ — 4 

F^-^T^ig E {lndetr(q,^.)-5]ln(4 6(«^) 



13=1 



+ (1- '^^o)ln 



(1 - e*"-'')(l - £-*"-'') 



}■ 



(54) 



where M has been taken, for convenience, to be an odd integer and the matrix T{q^LL}^) has elements [cf. Eqs. ( |45| ) 
and ^] 

r(9,^!^)a,a' = B{q,u^\a,a') +C{q,uj^\a,a') 

+B{-q, -ujy\a' , a) + C{-q, -uj^\a' , a). (55) 
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The first term within braces in Eq. ( p4[ ) has been obtained by the method discussed in Appendix ^ , which is required 
since the relevant Gaussian matrix is not Hermitian. The remaining terms within braces in Eq. (|54 ) originate instead 
from Jacobian contributions which are specific to the "radial" gauge and are included here for convenience of the 
following calculation. Note also that the frequency sum in Eq. ( |54|) has been symmetrized between positive and 
negative values, in order to extract directly its "continuum limit" (see below) according to a theorem reported in 
Appendix 

The mean-field contribution Fq to the (site) free energy given by Eq. (|5^) does not present any pathology in the 
i5 — > limit. One obtains the familiar expression for an (effective) system of independent particles 

BZ 

Jim^/?5^°^^--l5:5:in(l + e"^-(^)) (56) 

[cf. Eq. (^^], as shown in Appendix The fluctuation contribution Fi to the (site) free energy, on the other 
hand, has pathologic behavior in the S —* limit owing to the presence of the bosonic frequency sum in Eq. ([54|). 
It turns out, in fact, that performing the Matsubara frequency sum with the discretized time mesh (i. e., keeping 
< (M — l)/2 and taking the M ^ oo limit only after having evaluated the sum) gives a different result from the 
one obtained by interchanging the continuum limit with the frequency sum. The latter procedure requires standard 
mathematical techniques and has invariably been considered in the previous slave-boson literature. |— ■ 

The pathologic behavior of a bosonic frequency sum of the type (|5^) has been discussed in detail in RefM in the 
context of the simpler U = oo one-slave-boson problem. Specifically, it has been proved that this frequency sum can 
be suitably partitioned into a "continuum limit" contribution (where the (5 — > limit is taken at the outset in each 
term of the sum while the Matsubara index v is extended from —oo to +oo) plus a "contribution from infinity " 
which accounts for what is left out by the continuum limit . For the reader's convenience, we report in Appendix ^ 
the statement of a theorem which shows how the "contribution from infinity " can be extracted for a class of bosonic 
frequencv-sums that satisfy certain requirements. For a complete proof of the theorem we refer instead to Appendix 
B of Refill. p 

It has been also demonstrated in Ref.Efl in the context of the U ~ oo one-slave-boson problem that taking into 
account the "contribution from infinity " (to the bosonic frequency sum entering the Gaussian correction to the free 
energy) is not only required at a formal level, but it is also essential in practice for validating the 1/N expansion for 
the free energy. In the next Section we will show that taking properly into account the "contribution from infinity " 
to the sum (|5j) is essential for obtaining meaningful results at the Gaussian level even for the four-slave-boson case. 
This will hold in spite of the nontrivial additional complications related to the presence of the subsidiary operator R 
in Eq. (|. 

In essence, the reason why taking the continuum limit at the outset in the frequency sum of Eq. (|^) leads to 
incorrect results, is that the bosonic frequency uj^, enters the relevant fluctuation matrix (|5^) only through the phase 
factor exp{ia;i/(5}. By taking the continuum limit , one effectively regards the argument of this phase factor to be 
small (i. e., \uJi,\6 <C 1) and accordingly replaces exjp{iuj,^6} —f 1 + iuji,S. The flaw of this procedure is that, even in the 
limit S ^ 0, \uj,^\S can be comparable to unity since, by definition, oj^S = inv/M {v = —{M — l)/2, • • • , (M — l)/2) 
and max jw^Jj = tt. The replacement exp{zci;^(5} ^ 1 + iuj,j5 is thus not allowed to obtain the correct value of the 
bosonic frequency sum in Eq. (^4|) , and the limit 6 —> Q can be safely taken therein only after having performed the 
frequency sum. 

To be more specific, a bosonic frequency sum of the type (|5^ ) can be effectively partitioned into two contributions 
associated with "small" and "large" frequencies, respectively. The contribution from "small" frequencies (such that 
\uj^\6 ^ 1) corresponds to the usual continuum limit. The contribution from "large" frequencies (such that \ijJu\5 is of 
order unity) requires instead expanding the argument of the sum in terms of all 5 factors not entering the combination 
exp{ia;^(5}. Specifically, the contribution from "large" frequencies we are concerned with originates from terms 0{5) 
in this power expansion, since the sum of a large (of order M) number of terms each 0{5) yields a finite contribution. 

According to the general procedure of Appendix ^ extraction of the contribution from "lar ge" frequencies (that 
we have named "contribution from infi nity " ) amounts to isolating the function g of Eq. (|B10| ) and to determining 



its constant term 170 of the expansion (B9). Whenevet_the terms in the bosonic sum have simple expressions (as 
for the U = 00 one-slave-boson problem treated in Ref.tj), extraction of the "contribution from infinity " go can be 
done analytically with moderate effort. In the present four-slave-boson case, however, the expressions of the matrix 
elements (^5|) are exceedingly complicated to handle the sum in Eq. (jsj) analytically. For this reason, we have 
developed a computer algorithm for symbolic calculations which (i) evaluates explicitly the quantities (p3|)-(|35|) in 

terms of the bosonic mean-field values ag""* (a = 1, ■ • • , 8) and of the parameter S, (ii) sets up the matrix elements 
( p5| ) , and (iii) performs the steps indicated in Appendix ^ to extract the "contribution from infinity " g^. 

In this context, we found it convenient to regularize at the outset the matrix (|5^) for ^ in the limit of small S, 
by introducing an auxiliary fiuctuation matrix F such that 
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(57) 



where Xa — 5+i- This definition makes the matrix F and its inverse regular in S about 5 = and impHes 

that det r — S^detT. Accordingly, we write (in matrix notation) 

T{q, C\S) = To{q, + STiiq, C) + 0{S^) (58) 

with C = exjp{iuj^6}, and expand 

Indetr = tr Inf = In det To + Str (ro^^ri) + 0{S^) (59) 
provided Tq is nonsingular. For the lowest-order term we obtain 



detro(C) = (l-C)(l-r')n4&o''^'' 



(60) 



which (together with the factor from the transformation (^7|)) cancels the Jacobian contributions on the right-hand 
side of Eg. ([5^) (except for a term — ^ In 4 ^ that will be added to the u = Q component of the fluctuation matrix). 
We are thus left with extracting the constant term from the function 



5q(C) = tr [ro(C)-^ri(q,c)] 



(61) 



for any given wave vector q, according to the prescription discussed in Appendix 

We eventually obtain the following expression for the "contribution from infinity " to the Gaussian free energy (p4): 



dFo 



2+ 



BZ 



(62) 



which holds specifically for the paramagnetic case whereby Ag = q ^^'^ = £a{k) are independent of tr. In 
Eq. (^), fpis) is the ordinary Fermi function (which is recovered in the present context as the continuum limit of 
the function /m(£) - cf. Appendix P), Fq is given by Eq. (^, and J^lTZibg] depends on the specific form of the 
subsidiary operator R in Eq.(^). In particular, ^[7?.;6q] vanishes in the simplest case when R= 1. In all other cases, 
J-[TZ;bQ\ depends in a nontrivial way on the subsidiary function TZ introduced in Eq. (26) and on its first and second 
derivatives, which have to be calculated with the arguments (^^ taken at the mean-field level. Specifically, JF[7?.;&g] 
can be represented in the compact form: 



(1) ^(2) ran ^d^^ }_dn ^ ion 

\dni dn2 2 9713 2 dn^ 



(3) (4) /I 97^ 1 dTZ on dn 



In this expression, TZ ~ TZ{ni, ^2,^3, 714) with Ufj = &o'^''^ (/? = 1 



(3=1 

•,4), 



9^ 

dnl 



^0 = (4^'' ^0^^ + &o^^ foo'*')7?.(ni,?l2,n3,"4) 



(63) 



(64) 



is the mean-field value of the bosonic operator Zi^a given by Eq. (^, and b^^^^ = b)^' for the paramagnetic case we 
are considering. 

Equations ( p2| ) and (|6^) are the main results of this subsection. It will be shown in Appendix ^ how these results 
can be obtained alternatively within the "Cartesian" gauge. The reason to work here with the "radial" gauge is that 



u(3) 
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the resulting "continuum limit " calculation is free from infrared divergences. In particular, the "continuum limit 
contribution to the Gaussian free energy is given by 



1 



2/3Af 



BZ +00 

E E <|i- 



det Tc{q,uJi,) 



•5^0 + (1 - Cfi 



K 



-EM4 6r) 



/3=1 



(65) 



where rc(q,w^) is the continuum [5 0) limit of the matrix ( ^5| 

The physical free energy (per lattice site) at the Gaussian level is then given by [cf. Eq.(|5: 



Ad) 



(66) 



where F^'^'' and F^"'' are given by Eqs. 

by taking the mean-field values for the bosonic variables fop- In this way, each term dFo/dbQ^'^ in Eg. (|62|) vanishes 

The "contribution from infinity " F^"^^ in Eq. ( |66| ) is what has been missed in the previous literature on the 
four-slave-boson approach to the Hubbard model when considering Gaussian corrections beyond mean field. We have 
attributed the reason for this omission to an improper handling of the continuum limit of the functional integral. We 



Ad) 



) and (p3), respectively, with each contribution evaluated at self- consistency 



now argue that the occurrence of F^ 
within the functional integral. 



(d) 



is intrinsically related to the need of recovering the correct bosonic commutators 



D. INTERPRETATION OF THE "CONTRIBUTION FROM INFINITY" 

The "contribution from infinity " (|6^), when taken at self-consistency, can be interpreted in an euristic way by the 
following argument. Let's consider first the simplest case when 71=1 (and J^[7^;6q] = 0, according to Eq.(|6^)). In 
this case, noncommuting bosonic operators enter the Hamiltonian only through the number operators b^^"^^ lr\ We 
argue that taking naively the continuum limit of the action at the outset corresponds effectively to replacing each 
number operator in the original Hamiltonian as follows 

biP)Hm ^ i(J,(/3)t5(/3) + ^(/3)5(/3)t) = ^m^iP) + ^,(/3)t] ^ (67) 

and then treating the action associated with the modified Hamiltonian correctly (that is, by keeping the discretized 
time mesh until the end of the calculation) . To restore the correct original Hamiltonian, one has thus to subtract one 
half of the commutator whenever the term occurs in it. However, since bosonic commutators 

contribute terms of order at least 1/iV to the action of the functional integral (cf. Appendix in particular Eq. (A3) 



), one half of each commutator has to be subtracted only when fluctuation (1/A^) corrections are included. Thus, to 
the terms 

a 

-E ^lM,aS^,a + d\di) + U d\d, (68) 
a 

at a given lattice site in the original Hamiltonian, there corresponds subtraction of the terms 

Agi4-2Agi2 + C/i = 2AS-2A[j + ^ (69) 

from the Gaussian contribution to the (site) free energy calculated in the continuum limit (in the paramagnetic case) . 
In this way, we can account for the presence of the first three terms on the right-hand side of Eq. ( |62| ) . 

In the general case when TZ ^ 1 (and JF[7?.;5q] is nonvanishing) , the last term on the right-hand side of Eq. ( |6^ ) 
originates from additional noncommuting operators which are present in the bosonic (hopping) operator (^ as soon 
as the subsidiary operator i?^^^ is different from unity. In this case, noncommuting bosonic operators enter the 
Hamiltonian also through more general combinations than the number operators 

Quite generally, we argue that taking the continuum limit of the action at the outset corresponds to replacing the 
original Hamiltonian operator by a different operator Qc{H), obtained by a suitable average over the permutations 
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of the sequence of creation and destruction bosonic operators defining H. To restore the original operator _ff , one has 
then to introduce the difference Coo{H) between H and Qc{H) and set 

H = QciH)+Coo{H) . (70) 

By construction, the difference Cqo{H) is expressed in terms of bosonic commutators. It thus contributes to the action 
only terms of order at least In our case this implies that Coo{H) contributes to the Gaussian correction to the 

free energy only its average value C^{H) evaluated at the mean-field level. 

We can formulate the following hypotheses on the operator transformation Qc acting on a generic operator P: 

(i) Qc is a linear transformation, in the sense that Qc{P + P') = Qc{P) + Qc{P') for any given pair of operators 
P and P'; 

(ii) When two operators P and P' refer to "different kinds" of particles (i. e., to different bosons at the same or 
different sites, or to the same bosons at different sites), then Qc{PP') — Qc{P)Qc{P')', 

(iii) If K is a c- number term or a fermionic operator, then Qc{KP) ~ KQc{P)^^ 

With the above hypotheses, we can restrict ourselves to considering the transformation Qc acting on a monomial 
operator P entering the expression of the Hamiltonian operator i7, of the form 

P = (71) 

(with m and n arbitrary integers) for any given kind h of bosons. Accordingly, we might take Qc of the form 

Q,(5t"5™) - — y 7'fc(6^"6'") (72) 
n-p ^ — ' 
{fe} 

where k labels the n-p = [n ^ m)\ / {n\m\) permutations Vk of the product Alternatively, we could consider the 

simpler form 

Q^(6t»5m) ^ i(5tn5m + fe^fot"). (73) 



For our purposes, the two operators ( |72[ ) and ( |73| ) are equivalent to each other, since the average (mean-field) values 
C^(P) of the corresponding operators Coo(P) differ by terms of order at least {1/NY with respect to the average 
value of In the following, we shall thus consider only the simpler form (|73|). 

To evaluate the average (mean-field) value C^(6^"6™), we exploit the relation 



from which 



[6t,6™] = , (74) 



1 f),j.n am 



1! dx '^=1 dy ^y=^ 
2! dx^ '^=1 dy^ 'j'^i 
3! dx^ 1^=1 V l^'^i 

(75) 



Thus: 



where the remaining terms originating from the commutator (|7^) have been omitted since they contribute to (P) 
at order at least {l/N)^. At the relevant 1/N order, comparison of Eq. ( [76|) with Eqs.(^o|) (written for P in the place 
of H), (|7^), and ^) yields eventually 

Coo(&^"&™) = _ (77) 
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The corresponding average value is then obtained from Eq. ( [77| ) by replacing b and 6^ by their mean-field value 
With the monomial (^ij) we can construct the generic (normal-ordered) operators of interest, namely, 



: /(6tfe) 
: f{b^b)b 

The corresponding values of are , in the order: 

C^O fib^b) :)= 



J2 

oo 

J2 fnb^^^'b", 

oo 



(78) 
(79) 
(80) 



n ,2(«-i) 



1 



Cl(:6t/(6t6):)^„5o^/, 



(n + l)n 2(n-i) 

o 



bo 



{2f'ibi) + biribi)), 



Cl{:f{b^b)b:)^Cli: b^f{b^b):) 



(81) 



(82) 
(83) 



where we have assumed the series X^i^o Z"^" represent the analytic fmiction f{x). 

Thus far we have considered a single kind of bosons. For different kinds of bosons, can be evaluated by exploiting 
the following property. Given two operators Pa and Pb referring to bosons a and b, respectively, hypothesis (ii) above 
together with the definition (|70| ) imply the identity: 



CooiPaPb)= PaPb - Qc{Pa)Qc{Pb) 

= PaPb - {Pa - Coo{Pa)){Pb - Coo(Pb)) 

= PaCooiPb) + Coo{Pa)Pb - Coo{Pa)C^{Pb) ■ 

Taking the average value of both sides of Eq. (|8^), we obtain at the leading order in 1/A^ 

Ci^iPaPb) =< Pa >0 Cl,{Pb)+ < Pb >0 C°^{Pa) 

where < P >q signifies that the average value of P is evaluated at the mean-field level. 
The compound operators of interest are of the form 



: fia\a,b\b) := : a^''a''ia^ar{b^brbH^ 



(84) 
(85) 

(86) 



where the integers /i, v, p, and r can take the values and 1. In particular, when ii = v = p = t = wc obtain from 
Eqs. (H) and (IttI): 



(: /(ata, 6t&) :)= ^ (a^ (&^"&") + (a^-a™)) 



d d'^ 



dy dy^ 

f d d'^ \ ( d 9^ \ 



Cl{:f{a^a,bl):) + Cl{:f{al,b^b) :) 



(87) 
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where use has been made of Eq.(|l|) and of the assumption that the series J2m nfmnX"^y^ represents the analytic 
function f(x,y). By the same token, we obtain for the other cases of Eg. (pq) : 



f{a\a,b\b) :) = C^(: fiao,ao,b\b) :) + f{a\a,bo,bo) :) 



(88) 



Extension of Eq. ( pq ) to several kinds of bosons is straightforward. 

We are now in a position to reproduce the last term on the right-hand side of Eq.(p2[). This term turns out to 
be the contribution originating from the presence of noncommuting bosonic operators in the kinetic part of the 
Hamiltonian, which occurs only when nontrivial forms of the subsidiary operator in Eq.(^ arc considered (i.e., when 
Ri,<T 7^ !)• According to the above prescriptions, we then obtain for the kinetic part of the Hamiltonian: 



t 



E {fLn+^,^)^^o (C^(: zl :) + C^(: :)) 



(89) 



where Zq is given by Eq.(|64|). In deriving Eq. (^9|) we have made use of Eq. { p5\ ) and of the fact that, at the order 
of the 1/N expansion we are considering, the fermionic average has to be evaluated with the bosons taken at the 
mean-field level. Equation ( ^ ) further simplifies for the homogeneous and paramagnetic mean-field solution we are 
restricting to. In this case 



C^^0 4.0 = C^L(:z.+A,. ■■) = Cl{:z, :) 



(90) 



is independent of i and a. [It is here understood that the operator corresponds to a given reference site, say io]. 
We obtain eventually: 



2t 



(/i>/j+A,, 



i,A,( 
BZ 



'Af 

with 7(fc) given by Eq.(^l]). 

We are thus left with evaluating C^{: Z| :). Taking Z| of the form [cf. Eq. 



with 



by our conventions, and using Eqs. (^l|)-(p3|) and (|88|), we obtain 
where [cf. the dictionary (Esl)] 



c: 



'^(: 4^^^ :) = Ct {b^^' : 7^(6(^)t b^^\ b^^\ b^ , b^) : b^^^) 



+^^^(4^^7^(6«^6(2)t5(2),5(^)^4^)^):6(^) 

+CS.(: b^^'^^nibi'^\ b^^^\ b^'^H^'\ bi'^') : bl,'^ 





■ 5 




2 ° ° \ 


dm 


° dn\ dn2 










i9n4 ° dn'l_ 



92 



7^ 



7^ 



(91) 



(92) 
(93) 
(94) 



(95) 
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with the same notation used in Eq.(|63[). The other contribution C^{: e'^R-^s^ :) can be obtained from Eg . jos]) with 
the replacement 1^4 and 2^3 for the boson indices. Comparison of Eqs. ( |9^ ) and ( p5|) with Eq.(|63|) identifies 
eventually 

.F[7^;5o]=2zoC•^(:zT 0- (96) 

This result, in turn, implies that the last term on the right-hand side of Eq.(|6^) coincides with the contribution 
( ^ ) from the kinetic part of the Hamiltonian. 

The method that we have here provided to rationalize the "contribution from infinity " (^2|) might also serve as a 
practical prescription to avoid the burden of keeping explicitly the discretized nature of the functional integral until 
the end of the calculation.E3 This method may thus be useful when considering Hamiltonians different from and 
more general phases (like magnetic ones). However, extension of this procedure to physical quantities other than the 
free energy and its derivatives (for instance, to the electronic correlation functions) is not straightforward and requires 
a separate study. This problem will be touched upon in Appendix ^ in the framework of the "Cartesian" gauge. 



E. EXTENSION TO OPERATORS NOT EXPLICITLY IN NORMAL-ORDERED FORM 



Thus far we ha ve co nsidered subsidiary operators Ri^^ that are explicitly written in normal-ordered form. In 



fact, in subsection II C the normal-ordering prescription en tered the operative definition of the subsidiary function TZ 
[cf. the comment following Eq. (^7|)], while in subsection |II D| the normal-ordering prescription for all terms in the 



Hamiltonian was assumed throughout [cf. the main result (|9q) therein]. Nonetheless, it might be also relevant to 
consider operators, like Rf'^ given by Eq.(|5b|), which arc not explicitly written in normal-ordered form. For instance, 
it will result from the numerical calculations presentedpiii the next Section that the choice ( [5b| ) remedies the unphysical 
results obtained with its normal-ordered version (^^ O 

It is clear that, when adopting a functional-integral formalism based on coherent states, any non normal-ordered form 
for Ri^„ should be preliminarly rearranged into a sequence of normal-ordered terms. In general, this rearrangement 
might result into hardly manageable expressions (or even into divergent expressions when non polynomial operators 
like (^) are adopted) . This difficulty can be overcome as follows in the spirit of the 1 /N expansion considered in this 
paper. 

Let's consider first the rearrangement of the simple bosonic monomial operator (b^b)" into a sequence of normal- 
ordered terms: 

n 

{b%Y = Y,an,k ■ {b'^hy-^ : . (97) 

It is clear that a„_o = 1. A recursive relation can be readily obtained for the other coefficients of the expansion (p7|): 

a«,fc = a«-i,fe + (^^ - A:)a„_i_fc_i (98) 

with the initial conditions a„^„+i = a„._i — and ao,o = 1- lu particular, Eq. ( p8| ) gives for fc = 1 

nln — 1) , , 

an,i = ^ ^ ' ■ (99) 

Knowledge of coefficients with A; > 1 is not required at the order we are considering in this paper. 
For a general bosonic operator we then write 

^ ^ n=0 ^ ^ n=0 fe=0 

where the expansion ( |97| ) has been used and b^b has bee n div ided by N according to a standard procedure of the 1/A^ 
expansion [see Appendix M . The right-hand side of Eq. ( 100 ) admits a direct representation in terms of the functional 



integral, by mapping : {b*bY : onto {b^bm-iY in the action for any (positive) integer q. In the spirit of the 1/7V 



expansion, we then rescale the bosonic variables bm by setting bm = \Nbm- The operator (100) is thus mapped in 
the action onto the series X]fc!Lo ^"'^Z''^'' (^m^m-i)' where 

oo 
q=0 
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At the first order in 1/iV, only terms /^'^•'(a;) with fc = and 1 have to be retained in the a ction. Specifically, the term 
with k = will c ontribute both to the leading term Fq of the free energy [cf. Eq. ( |A10| )] and to its correction 
Fi [cf. Eq. (A13)], by considering its mean-field value and its Gaussian terms, respectively. The term with fc = 1 will 
instead contribute only to -Fi, by considering its mean-field value.cj In this way, the contribution to Fi of the term 
with fc = 1 can be effectively reabsorbed into the "contribution from infinity" to Fi originating from the term with 
fc = 0. Specifically, it turns out that this contribution simply cancels the term containing the second derivative of / 

in Eq. i^. 

Entering Eq. (M) into Eq. (101) we, in fact, obtain for f^^\x): 



9=0 



1 



xf'ix) 



SO that 



(102) 



(103) 



owing to Eq. (81). In Eq. (103) we have consistently replaced the argument of the function f^^^x) by 6q and set 
= 1 eventually. Generalization of the above results to the case of several bosons is straightforward. 

In conclusion, for an operator not explicitly in normal-ordered form the above considerations have the effect of 
modifying the factor JF[7?.; 6g] in Eq. (|62|), by (i) dropping the last terms in Eq. (63) containing the second derivatives 
of TZ and (ii) interpreting the function TZ in the remaining terms in Eq. ( |63[) (as well as in the expression for the 
continuum part of the fluctuations) as being associated with : Ri „■ :, i.e., with the normal-ordered component of Ri,^- 

It will be shown in the next Section that, when specified to R^^^ , the above prescription succeeds in eliminating the 
unpleasant features occurring for its normal-ordered version Rf^^ [cf. Eq. (|5a|)] . Namely, it eliminates an unphysical 
divergence of the ground-state energy when the zero-occupancy limit is approached (making the ground-state energy to 
correctly vanish in this limit) and completely suppresses the fluctuation corrections to the KR mean-field ground-state 
energy for U = and any filling value. 

Finally, we remark t hat the right-hand side of Eq. ( |103| ) could be obtained directly by generalizing the euristic 
argument of subsection |IID| to the case of non normal-ordered operators. To this end, one suitably defines the operator 
Qc{P) (and consequently Coo{P)) for a non normal-ordered operator P by extending the "specular" prescription (|7^), 
and evaluates Cj^ accordingly. By this procedure, one is led to interpret the right-hand side of Eq. (103) as the 
"contribution from infinity " originating from the non normal-ordered operator f{b^b). 



III. NUMERICAL RESULTS FOR A ONE-LEVEL TWO-SITE MODEL 



In the previous Section we have shown that the additional contributions (62) to the Gaussian free energy originate 
from a careful handling of the continuum time limit of the functional integral. The question naturally arises whether 
these additional contributions might substantially affect the numerical value of the free energy (and of related physical 
quantities) in the practical cases of interest, or even modify its behavior in a qualitative way. Specifically, we may 
inquire whether the additional contributions (^2|) could serve to overcome the inconsistencies encountered when 
including Gaussian fluctuations by treating the functional integral in the continuum limit . In this respect, we have 
found that the wrong curvature of the free energy versus U results with the conventional continuum limit treatment, 
thus producing an unphysical increase of the (average) number of doubly occupied_sites for in creasing U. Additional 
inconsistencies have been pointed out for the noninteracting ([/ = 0) case in Ref.E3. 

In this Section we will show that inclusion of the "contribution from infinity " ( p2[ ) indeed succeeds in overcoming 
the inconsistencies mentioned above. In addition, we will show that relaxing the normal-ordering prescription ( |5a| ) of 
the KR choice for the subsidiary operator, makes the 1/A results coinciding with the exact (free-particle) solution at 
U = for any filling and number of lattice sites. In this way, the goodness of the KR mean- field solution will not be 
spoiled by the inclusion of fluctuations. 

We shall illustrate the above effects by numerical calculations at a generic value of U for a one-level two-site 
model, which avoids the complications due to the spatial structure while dealing explicitly with the imaginary time 
discretization. We expect, however, these effects to hold for a more general multi-site system as well, since the 
time continuum limit does not depend on the spatial structure [In this respect, we will also present some numerical 
results for the lattice case when U = 0]. The restriction to a two-site model will also enable us to compare the 
numerical results, obtained withip-altcrnative approximations to the functional integral , with the available exact 
solution discussed in Appendix In addition, to get a practical insight into the usefulness of a given choice for the 
bosonic hopping operator z, we compare the results obtained with several alternative choices of z. 
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A. MEAN-FIELD SOLUTION 



At the mean- field level , the free energy (per lattice site) of the paramagnetic solution is given by [cf. Eqs. (|5^), 
@, and (H)]: 

Po = --,j2^n{l + e-^^^)+Fl,^^ (104) 



' p=i 



where 



and [cf. Eq. (M)] 



) = XUb'i^' + 26(^)2 + - 1) - 2Ag(6(2'^ + b^^^') + Ub^^^' (105) 



£1 = Af; - - 2tpl 
£2 = A[5 - /i + 2tpl 



(106) 



with the notation 



PO - b''o\b[P + 5[,^))7^(71l, 712,^2,714) (107) 



in the place of zq [cf. Eq. (|6^)]. We have here specified to the two-site model the general expressions reported in 
Section || for the lattice model, and made use of the paramagnetic ansatz b^^"^ — 6g'^^ . 

The parameters b^\ bl^\ b'^\ Aq, Ag, and p are determined by minimizing Fq + fin (where n stands for the particle 
density per lattice site). The resulting mean-field equations are: 



4t(/F(£2)-/F(£i))po^ + ^^=0 (/3 = 1,2,4), (108) 



db)^> dbl"" db'o 

^=6r + 26r+fe^'''-l = 0, (109) 
dFp 



/f(£i) + /F(e2) - 2{bl^^' + 6(^)2) = , (110) 



r) F 

= friei) + fFie2) = n . (Ill) 



In particular, in the zero-temperature limit Eq.( lll ) implies 

/F(ei) = l , /F(£2) = n-1, (112) 

whenever n > 1 and po is nonvanishing, with the chemical potential jumping from ei to £2 across n = 1. 

At "half filling" (i. e., when n — 1) there exists a critical value Uc such that 6g^'' — b^^^ = and fog^' = l/\/2 for 
U > Uc (and the associat ed H cl mho ltz free energy density. Fq + pn vanishes in the zero-temperature limit). It can be 
readily shown from Eqs. ( |108D -( pTl| ) that Uc is given byEJ 

c/c = 4^7^2(o,i i 0) . (113) 



In deriving Eq.(113) we have exploited the symmetry properties 

7?.(7il,7l2,n2,«4) = 7^(714, 7l2,7l2, Til) = 7^(7l2 , Tll , 714, 77.2) (114) 

whi ch ar e satisfied by the choices we shall consider for TZ. This transition a.t U = Uc makes the effective hopping in 
Eq.( |l06|) vanishing, and it is referred to as a Mott-Hubbard transition.! For the two-site model we are considering 
this transition is clearly an artifact of the slave-boson approach at the mean-field level , which will be (at least partly) 
cured by the inclusion of (Gaussian) fluctuations. 

For comparison, we also give the results of the conventional Hartree-Fock decoupling. For the one-level two-site 
model of interest, the Hartree-Fock ground-state energy (per lattice site) in the paramagnetic phase is given by 
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E"^ = 2t{f,{eD fF{er)) + (115) 



where 



-HF ,, 



with fpi^i) and /i^(e2) given by Eq. ( |112| ) (in the zero-temperature hmit for n > 1). Since e^^ (p = 1,2) are 



i ndep endent of U (at fixed particle density), the energy ( 115 ) grows hnearly without bound for increasing U. Equation 



(115) also shows that, within the (paramagnetic) Hartree-Fock decoupling, the site probability of double occupancy 
is given by (n/2)^, irrespective of U. One thus concludes that the (paramagnetic) Hartree-Fock decoupling can be 
approximately correct only for U t (being, in particular, exact when J7 = 0). Numerical comparison of the mean- 
field ground-state (site) energy E = Fq + fin (in the zero-temperature limit) and of the associated 1/A^ fluctuation 
results with the Hartree-Fock E^^ and with the exact solution will be presented in the following. 



B. OVERCOMING THE PROBLEMS ASSOCIATED WITH THE CONTINUUM FLUCTUATIONS BY 

THE DISCRETIZED FLUCTUATIONS 



A naive handling of the functional integral (whereby the continuum limit is taken at the outset) yields, at the 
Gaussian level, the (site) free energy (|6^) with the term F^'^'' omitted and F^'^'' given by Eq.(|65|). In that expression, 
the continuum fluctuation matrix Tc{q, lOv) is obtained by performing the 6 ~* Q limit of the matrix (|55| ) (with B and 
C given by Eqs.(^5|) and (^), in the order), i. e., by letting (5 — > everywhere 5 appears. 

In particular, the continuum limit of the matrix B given by Eq.(^) has elements 



Bc{q,uj^\l, 1) = iuj^ + Ag, 

6,(9, c^, 1 1,2) = ••• =6,(q,c^,|l,7) = 0, 

Sc(9,c^.|l,8) =6c(g,t^.|8,l) 



(117) 



and so on; the continuum limit of the factors S^^'^ entering the expression ( |5l| ) for the matrix C is given by 

£(i)(g,c.,;fe,r |/3)=teo(7(fc) +7(^-9))^, (/?=l,---,4) 
£(1) (q, c^,; fe, t |5) - tz^n{ -f{k) - j{k 



(118) 



and so on (with zq given by Eq. (|6j)); the continuum limit of the factors entering the expression ( ^l|) for the 
matrix C is given by 



2.o^7(fc) 



dzo 



(1)2 



(1) 



(7(fc - q) + -f{k + q) ' 



(119) 



and so on.t3 Finally, the continuum limit of the Fermi function fM{e) and of the polarization function nM(e,e'; v) 
entering the expression of the matrix C are given in Appendix ^ 

The above expressions hold for a general lattice model. The two-site model that, we consider for numerical calcula- 
tions is recovered by restricting the wave vectors to the two values and 7r/|A|.E2l 

Proper handling of the functional integral (with the continuum limit taken only at the end of the calculation) adds 



to the continuum limit (site) free energy Fo -I- f\^^ the "contribution from infinity " F^"' given by Eq. ( |62|) with 
the appropriate form (p3) for T\K.\hf^, with or without the terms containing the second derivatives depending on 



(d) 



the ordering prescription (cf. subsection HE). Our purpose here is to demonstrate to what extent the contribution 



f['^^ modifies the free energy in a quantitative and even in a qualitative way, using the two-site model as a simple 

prototype system. In this way, the contribution F^'^'^ will prove to be important also from a practical point of view. 

To proceed further, we have to select an explicit form for the subsidiary operator Ri „ introduced in Eq.(Q). The 
choice commonly adopted in the four-slave-boson literature is the original KR choice (^). The associated subsidiary 
function introduced in Eq.(pq) takes the form 
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7^^^(ne, n,, ns, nj,) = ,^ (120) 
VI - "e - riffVl - rid - no- 

which recovers at the meaiL-field level the (paramagnetic) Gutzwiller approximation for the lattice Hubbard Hamil- 
tonian at any band fiUingo However, it will turn out that the normal-ordering prescription leads to results which 
markedly depart from the free-particle solution at J7 = even when the "contribution from infinity " is properly in- 
cluded. In addition, this prescription leads to a divergent free energy in the zero-occupancy limit. These shortcomings 
can be remedied by suitably relaxing the normal-ordering prescription, as shown in the next subsection. 

Figure |l| compares the ground-state energy per lattice site (in units of t) for the one- level two-site model versus U/t 
at "half filling" (n = 1.0) and at ti = 1.2, as obtained from the exact solution of Appendix |^ and from alternative 
(param agn etic) approximations to the four-slave-boson method with the KR form (^^ for Rijj£3 The Hartree-Fock 
result ( |115| ) is also shown for comparison (HF). The KR mean-field solution (KR) is seen to be in good agreement 
with the exact solution (EX) for both doping values, while the fluctuation results considerably worsen this agreement 
regardless they are obtained by the incorrect continuum (CFL) or by the correct discretized (DFL) time limiting 

procedure in the functional integral (i. e., without and with the inclusion of the term F^'''^ of Eq. (|6^)). The HF 
result, on the other hand, is in agreement with the exact solution only for small values of U (say, when U ^ 5t). 
Note from Fig.|l|a that at "half filling" both the CFL and DFL curves terminate at a critical value of J7 ( ^ 6.8t) 
where an antiferromagnetic instability develops at the mean-field level , thus preventing the inclusion of paramagnetic 
fluctuations. Note also the arrow which locates the critical value Uc/t corresponding to the Mott-Hubbard transition 



for the KR mean-field solution (the choice ( |120| ) for TZ gives Uc/t = 16 according to Eq.(113)). 

We conclude from Fig.|l| that, with the KR choice ( pa] ) for Ri^a-, the DFL results can be even worse than the CFL 
results, as they depart from the exact solution more than the mean-field results. This finding could appear surprising, 
since one would expect a more complete treatment of the functional integral to produce better results. That this is not 
the case with the KR choice (js^) can definitely be concluded by looking at the DFL free energy in the zero-occupancy 
limit, which diverges for this choice in this limit. 

Before searching for alternative prescriptions for the subsidiary operator Ri^^r in Eq.(Q), it is relevant to verify 
whether the above unpleasant results were peculiar to the finite-size model. To this end, we have performed the 
U ^ calculations also for the lattice case versus the doping value n — 1 (since in the lattice case the antiferromagne- 
tic instability occurs already at infinitesimal U due to the perfect nesting of the Fermi surface). The results shown 
in Fig.^ confirm our conclusion that, regardless of its successes at the mean-field level , the KR form of Ri^a is not 
suited for the inclusion of fiuctuations. Note, in particular, the divergence of the DFL results when n — 2, which was 
mentioned above for the two-site model. 

This conclusion calls for alternative expressions of the subsidiary operator. One could initially resort to the simplest 
possible expression (||) for Zi^, for which Riu — R^ „ — 1 and TZ = 1. In this case the factor JF[7?.;6q] of Eq.(|6^) 
vanishes identically [cf. Eq.(|63|)]. The corresponding results for the ground-state energy (per lattice site) of the one- 
level two-site model are shown in Fig.||, with the same notations of Figjl] (except for the labeling of the mean-field 
curve, which is now MF). Note that the critical value for the Mott-Hubbard transition to occur at "half filling" is 



now Uc = 4t, according to Eq.(113). We see from this figure that the DFL results now improve the agreement with 
the exact solution with respect to the MF results, while the CFL results considerably worsen the agreement. This 
finding thus gives us support for treating correctly the time continuum limit of the functional integral within the 
four-slave-boson approach. 

In addition, we note the following shortcomings of the CFL results: 

(i) At finite doping (n > 1) the CFL result yield the wrong curvature of E/t versus UJi^ thus producing an 
unphysical increase of the (average) number of doubly occupied sites for increasing U /i.Lj This shortcoming is 
remedied by the DFL results. 

(ii) As mentioned above, in the zero-occupancy (n = 0) limit, the ground-state energy has to vanish irrespective of 
the value of U . The CFL results are instead finite in that limit. (This shortcoming is also shared by the previous 
KR choice for Ri^a). The mean-field results MF and HF behave instead correctly in the limit. Once more, the 
DFL calculation remedies the shortcoming (with the exception of the previous KR choice for Ri ^r which leads 
to divergent DFL results in the n = limit). The n — Q limit thus provides us with a common reference level 
for comparing alternative results for the ground-state energy. 

Although the choice R^ ^ — \ cures formally all fiaws occurring with the KR choice when fluctuation corrections 
are included, it is seen from Fig.^ that the mean-field solution associated with i??^ = 1 constitutes a rather poor 
starting point for the inclusion of fluctuations. In this respect, even the HF mean field is closer to the exact solution 
for U/t < 8. In particular, the choice R^^ = 1 does not reproduce at the mean-field level the independent-particle 
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result when L'^ = 0. By contrast, the KR form for Ri a- has been tuned to reproduce exactly the L'^ = result at the 
mean-field level , by introducing the nontrivial factor TZ^^ just as a "normalization" factor. 

It is thus clear from the above discussion that a more sophisticated prescription than = 1 is required for the 
subsidiary operator, which would preserve the nice features associated with RfJ^ at the mean-field level and, at the 
same time, would avoid the unphysical behavior resulting from Rf^ when fluctuations are included. 



C. ALTERNATIVE SELECTIONS OF THE SUBSIDIARY OPERATOR R 



An appealing feature of the KR (paramagnetic) mean-field solution is that (besides recovering the U = OJimit) 
it reproduces the Gutzwiller approximation for the single-band (lattice) Hubbard Hamiltonian (at any filling) .□ This 
approximation is, in turn, expected to-capture the essential electronic correlations for the Hubbard model by reducing 
the site double occupancy at large UE3 One would then like to select the form of the operator Ri,cr in such a way that 
the Gutzwiller approximation is still recovered at the mean-field level (at least for some particular filling). Recall in 
this context that the subsidiary operator Ri ^ of Eq.(^ has to be equivalent to the unit operator in the relevant Fock 
subspace. 

An alternative selection for the operator Ri ,^ might rest on the following criterion which exploits the mean-field 



constraints (109)-( |lllD in the zero-temperature limit. In particular, at half filling [n — 1) the relation 6q + 

^o'^^^ + ^0^^^ — 1/2 holds for any value of U, so that any function TZ of 6g^^^ -I- 65^^^ and 6q^^^ -f 5q^^^ will also be 
independent of U . This property is satisfied by the KR form ( |12C| ). For any other form of TZ satisfying this property, 

it is then enough to "normalize" TZ in such a way that it coincides with TZ^^ when 6q^^^ -I- fop^''^ — b'^^^ + b'^^^ — 1/2. 
This criterion guarantees that this new TZ and TZ^^ give the same mean-field results. Note that the mean-field free 



energy Fq is also independent of TZ (whenever the property (114) holds) since Ag does not enter Fq at self-consistency. 

A form of the function TZ, which reproduces the_KR mean-field results at half-filling, is obtained by suitably 
"linearizing" the expression (|120|), namely, by settingE3 



K^^^ {ne,na,ng,nd) = [1 + x{ne + ng)][l + x{na + Ud)] (121) 
with X — Xq such that (in the paramagnetic case) 

This gives xq = 2(\/2 — 1) = 0.828 for the positive solution. Although this particular value of x has been selected 
at "half filling" , in the following we shall use the same value of x for any filling. Note that the subsidiary function 



(121) corresponds to the "linearized" subsidiary operator i?,f^^ = [1 + x{e\ei + s\ _^Si^-a)][\ -I- x{d\di -\- s\ ^Si^a)] for 
which_the normal ordering is irrelevant. This simple choice is free from the nonanalyticity problems affecting instead 



The results for the ground-state energy (per lattice site) of the two-site model corresponding to (121) are shown 
in Fig.|| for filling values n — 1.0, n = 1.2, and n = 1.5. Curves are labeled with the conventions used in Fig||. 
Note that the DFL results at "half filling" are now quite close to the exact solution for C//t < 8. In particular, the 
DFL fluctuation corrections do not change appreciably (i. e., within less than 1%) the exact solution aX U = 0. 
(We have verified that this finding remains true also for an infinite system). The agreement between DFL and EX, 
however, is not so good for U close to the critical value Uc (= 16i) for the Mott-Hubbard transition to occur in the 
mean- field solution. Note also the presence of a "cusp" about Uc in the fluctuation results. This "cusp" is of no special 
concern since it is a characteristic feature of the Gaussian corrections near the transition point where higher-order 
fluctuations are expected to be important. This "cusp", in fact, quickly disappears as soon as one moves away from 
half filling (cf. Figs.^ and|^c). Note finally that the exact solution for the ground-state energy approaches a linear 
trend with respect to U for increasing doping. As a consequence, the HF approximation becomes progressively more 
reliable for increasing doping. 

The behavior of the ground-state energy (per lattice site) versus doping n — 1 is shown in Fig. ^ for three character- 
istic values of U /t with the same "linearized" TZ^^^ used for Fig. ^. The CFL results, however, have been excluded 
from Fig. ^ as they depart markedly from the other results. Even in this case, the correct fluctuation contribution 
(DFL) give quite good overall results. Note that the piecewise linear behavior of the exact solution versus doping is 
due to the finite size of the system (as well as to the zero-temperature limit) , since the exact solution for this system 
is constructed by interpolating between the solutions with n — 1.0 and n = 1.5 and with n — 1.5 and n = 2.0. Note, 
however, from Fig.^ that the free-particle results a.t U — are reproduced by the mean-field solution and by the 
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fluctuation corrections at "half filling" only (and when n = 0.0 or 2.0), while departures from the correct results 
become more pronounced at finite doping. (Recall that the KR mean field provides instead - by construction - the 
correct U — results at any doping and thus the correct value for the compressibility). 

The above shortcoming is obvi ousl y due to the fact that the parameter x of the "linearized" TZ^^^ has been selected 
at "half filling" according to Eq.(122). Although one might fix x at any selected doping, it is certainly not satisfactory 
to adjust X to be doping dependent in order to reproduce the independent-particle results at any doping. For this 
reason the "linearized" form (121) proves to be not completely satisfactory. 

The fact that the KR mean field reproduces the {U = 0) free-particle results for any doping is due to the perfect 

in Eq.(107) with the self-consistent doping 



cancellation of the doping dependence of the factor ^q^' ( feQ^-* 



dependence of the TZ^^ factor (12C). It is also evident that no "polynomial" generalization of (121) can yield this 
perfect cancellation for all doping values. 

To fully preserve the nice features of the KR mean field, we finally consider the "square root" form Rf^ ( pb| ) 
for the subsidiary operator, whereby the KR normal-ordering prescription in ( ^ ) has been removed. Following the 
discussion of subsection II E , this can readily be achieved within the 1 /N expansion by dropping the last terms in 
Eq.(|63|) containing the second derivatives and using TZ^^ in the rest of the calculation. 

The results obtained for the ground-state energy (per lattice site) of the two-site model with the choice a — Rf^ 
are shown in Fig.|| versus U for n = 1.0 and n — 1.2, where the mean-field (KR) and the complete 1/A^ (SQ) 
calculations are compared with the exact solution (like in Fig.^, an antiferromagnetic instability develops at "half 
filling" in SQ when U = 6.8t). The corresponding results versus doping n — 1 are shown in Figj^ for three characteristic 
values of U. Note, in particular, from Figj^a that the 1/N fluctuation corrections are now completely suppressed at 
U = for any doping. We are thus able to recover exactly the U = solution not only at the mean-field level but 
also with the inclusion of fluctuations, without having to adjust any parameter (like, for instance, the x parameter in 
Eq. (121)). In addition, this finding is not limited to the two-site model but applies also to the lattice case, as shown 
in Fig.^ In this way, the non normal-ordering SQ prescription (5b) overcomes the difficulties signaled in Ref.ta, 
without the need of changing the functional form of the bosonic operator z at the leading 1/N order beyond mean 
field. Note, finally, from Figs.^ and |^ that the SQ prescription gives good agreement with the exact solution for the 
two-site model at any doping even at finite U. In particular, for the filled band (n = 2) the SQ prescription recovers 
the exact result E = U ior any U [or = for the empty band {n — 0)], in contrast with the KR normal-ordering 
prescription that gives divergent results as n 2 (or n — > 0). This divergence resides in the second derivative terms 
of Eq. (63), which are appropriately eliminated by the non normal-ordering SQ prescription. 



IV. CONCLUDING REMARKS 



In this paper we have demonstrated the importance of taking correctly the continuum limit of the functional 
integral for the four-slave-boson method by Kotliar and Ruckenstein, by deriving in detail the additional (ultraviolet) 
contributions to the free energy associated with a proper handling of the vanishing of the imaginary time step. 
Althaugh these additional contributions to the free energy have been already derived for the one-slave-boson case in 
Ref.til, the presence of several bosons and of a nontrivial bosonic hopping factor z (which is peculiar to the KR method) 
results into additional noncommuting terms in the slave-boson Hamiltonian, which, in turn, make the proper handling 
of the continuum limit of the functional integral quite more involved. Nonetheless, the rather elaborate mathematical 
derivations that resulted proved necessary to implement correctly the KR method and make it working properly in 
practice. We have also provided an "effective" rule to obtain the additional terms in the free energy directly from 
the original Hamiltonian, without the need of going explicitly through the limiting discretization procedure of the 
functional integral. The fact that this careful procedure has anyhow produced unphysical results at the 1/A^ order 
beyond mean field when adopting the conventional KR choice for the bosonic hopping factor z, has further required 
us to search for alternative forms (or prescriptions) for z. Specifically, we have shown that removing the KR normal- 
ordering prescription for z enables us to overcome the above difficulties. In particular, we have verified numerically 
that the non normal-ordering prescription succeeds in suppressing completely the 1/N fluctuation corrections to the 
ground-state energy at C/ = and for any band filling. We regard this finding to be rather compelling to remove the 
doubts which have been raised on jthe capability of the KR four-slave-boson method to yield meaningful fiuctuation 
corrections (see, in particular, Ref.E3). 

In this context, some additional features still need, however, to be verified. For instance, it should be verified 
that removing the normal-ordering prescription from the KR choice of z suppresses the 1/A^ corrections at [/ = 
also for dynamical quantities (like the fermionic Green's functions). In this respect, we have developed preliminary 
arguments which show that the non normal-ordered subsidiary operator (5b) reproduces the exact U = Q results (at 



24 



zero temperature) for any value-of N , thus suppressing the fluctuation corrections to the free energy and to dynamical 
quantities exactly when U = OEB 

Although the formal results obtained in this paper are valid even at finite temperature, numerical results have been 
restricted to the zero-temperature limit throughout. Extension to finite temperature is then in order. Besides, one 
may also consider multi-band Hubbard-type Hamiltonians (like the Emery model for a Cu02 layer), for which using 
our "effective" rule to obtain the additional terms in the free energy should be straightforward. 

The numerical results presented in this paper have primarily concerned the two-site model, although results have 
also been presented for the infinite lattice when U = Q. Restriction to t/ = is connected with the paramagnetic 
ansatz we have considered throughout, since the matrix of Gaussian fluctuations develops instabilities for J7 = 0^ 
and n = 1. Consideration of magnetic phases is thus necessary, at least near "half filling" . We mention in this context 
that the KR four-slave-boson method, which is not manifestly spin-rotation- invariant, might not properly account for 
transverse magnetic fluctuations at the order 1/N. To this end, it would be necessary to include additional "angular" 
variables associated with transverse fluctuations. This could be achieved by the spin-rotation-inwariant method of 
Rcf.Q where two additional "s" bosons are introduced at the outset. Resorting to the method of Ref.Q, in turn, requires 
one to assess preliminarly the effects of the proper handling of the time continuum limit of the functional integral at 
the order 1/N that have been discussed in the present paper. Work along these lines is in progress. 
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APPENDIX A: 1/N EXPANSION WITH THE FOUR-SLAVE-BOSON METHOD 

In this Appendix we show how a suitable expansion can be constructed for the four-slave-boson problem. The 
need to organize the theory by means of a expansion originates from the lack of an intrinsic small parameter, 
in powers of which a more conventional perturbation theory could be defined. In addition, the introduction of the 
parameter 1 /N enables one to establish a meaningful comparison of the results obtained via different representations 
of the integration variables in the functional integral (i.e., by using either the "radial" or the "Cartesian" gauges). 

As for the one-slave-boson problcm,til also in the present context the integer 2N represents the number of spin 
components of the fermionic field coupled to the slave bosons. Also in this case, it turns out that the mean- field solution 
becomes exact in the N ^ oo limit and that the leading corrections to the free energy are associated with the 
Gaussian fluctuations. For physical quantities different from the free energy, however, the leading corrections 
require also the inclusion of higher-order fluctuations (which is equivalent to a l/A?^ shift of the mean-field parameters). 

Since the physical case of interest corresponds to the value = 1, one may wonder whether including only the 
corrections to the mean-field solution might be s uffic ient to obtain in practice sensible results when A'^ = 1. In this 



respect, the numerical results presented in Section III for a simple model system give us confidence that extrapolating 
the 1/N results down to A^ = 1 caxi_,work rather accurately for practical purposes (provided the form of the subsidiary 
operator is chosen appropriately) 

A meaningful 1 /N expansion for the four-slave-boson problem can be introduced in the following wayiEll 

(i) To begin with, the fermionic spin label (which in the original action (^)-(|l^) was restricted to the two values 
±1) is formally extended to run over the 2N values S — an with a = ±1 and n = 1,2, ■ ■ ■ , N. 

(ii) In order not to proliferate the number of slave bosons in the action (|7|)-(p^) (thus keeping consistently the 
number of mean-field equations independent of N), we retain only two "magnetic" bosons Sa-=+i and So-=-i, 
where now Sa=+i is associated with fermions with positive {S > 0) spin projection and Sg-^-i with fermions 
with negative {S < 0) spin projection. This corresponds to generalizing the constraints (|3|) as follows: 



CT = ±1 

E flsf^^s - 4,aS^.. ~ 4d, = (a = ±1) . (A2) 



In this way, the "spin" label a of the bosonic variables and Zi^cr in Eqs. (||)-(|T^) is maintained. Note that 



the replacement 1 ^ A on the right-hand-side of Eq. (Al) is necessary since the mean value of the fermionic 
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term in Eq. ( |A^ ) is of order N. By the same token, the unity in the square roots in Eq. 
in its "Hnearized" version discussed in Section III) has also to be replaced by N. 



(and consequently 



(iii) To keep the mean-field bosonic parameters independent of N, we rescale the bosonic fields by letting 

^i.m. ^ ^ N Ci. 777, , 



(A3) 



and similarly for Si_a-,m (f — il) and di_rn- 

(iv) Requiring the kinetic term ( [lO| ) to be homogeneous in N of the same order as the other two terms (||) and (||) 
of the action, leads to suitably rescaling the hopping integral t depending on the choice of the subsidiary factor 
Ri^cr in Eq. 1^. In particular, t tjN^ when Ri^^ ^ I; t ^ t with the KR choice (|l|); t t/N^ with the 
"linearized" choice corresponding to (121). 



With the above prescriptions, the action Sm is generalized as follows: 



M-1 N 

EE 



1 cr^itl <i,j> 



I i,(Tn,m^t,m;j,rrL 



j .(7n,7n—l 



NB[b] 



(A4) 



where < i,j > limits the lattice sum to nearest-neighbor sites. G and B in Eq. (A4) are functionals of the bosonic 
fields b (including now the Lagrange multipliers) which, by our assumptions, are independent of N. [For the non 
normal-ordering prescription for Ri^a, however, Q is itself written as a power series in 1/A^ - see subsection [IE and 
footnote |2^. When = 1 the action (A4) consistently reduces to the form (|7|)-(p^) considered in the text. 

It is important to emphasize that the A^-component generalization of the slave-boson Hamiltonian (which has led 
to the action ( |A4| )) is equivalent to the A^-component generalization of the original fermionic Hamiltonian (Q) only 
for N = I. The resulting 1/N expansion acquires thus physical meaning only when A^ is eventually set equal to 
1. However, this should not be considered to be an inconvenience of the present approach since the A^-component 
slave-boson Hamiltonian has not been conceived to study a physical system with truly large spin; rather, the resulting 
1 /N expansion must be considered as a purely formal tool to control the approximations systematically. 

The lack of correspondence of the two ( pure ly fermionic and fcrmionic-bosonic) Hamiltonian problems when A^ > 1 
results from (i) the constraints ( |Al[ ) and (A2) no longer providing (in general) a one-to-one correspondence between 
the Fock spaces on which the two Hamiltonians respectively act, and (ii) the matrix elements of a given operator 
being different even for pairs of corresponding states. As an example, let's consider the case N = 2. In this case. 



to the original configuration /J+3/2/i^-i/2l^ instance, the constraints (Al) and (A2) associate the two distinct 



configurations (ije| /J _|^3^2/J-i/2l'-' ^^"^ '^1 -i-i'^l -i//+3/2'^i^-i/2l^ On the other hand, to the original configuration 
fj ^3/2fi +i/2\^ there corresponds the only configuration (2)~^^^(s| _|_j)^/J_(^3^2/J+i/2l^ -^^ [^^ote that the mapping 
between the two Fock spaces (i.e., with and without slave bosons) preserves the fermionic configurations (specified 
by the operators and /t, in the order) and adjusts the bosonic configurations to satisfy the constraints]. For the 
latter state the average value of the interaction term Udldi vanishes, while for the corresponding original configuration 
the average value of the original interaction term UJ2s>S' fi sfi S'f^^S' f^^s equals U (with S and S' taking the 2A^ 
values — A^, 1,+1,---, -I- A'). Obviously, when A^ = 1 both the one-to-one correspondence between states and the 

equality of the matrix elements is fully preserved. 

Returning to the action (A4), we carry ou t th e associated 1/A^ expansion in the usual way by first integrating out 
the fermion (Grassmann) fields. The action (A4) is then replaced by 



S 



(AT) 



-N{tr logGik; cr = +1] + tr logg[b; a = -1]) + NB[b] = NS[b] 



(A5) 



where the trace is taken over the indices (i, j, to) of Eq. (A4) and S[b\ is independent of A^. Note that 5*^^^ depends 
on N only through the explicit factor N on the right-hand side of Eq. (A5), as a consequence of the assu mptions 
(i)-(iv) spelled out above. [Apart from the non normal-ordering prescription for Ri^a, for which S[b\ of Eq. ([A^ ) is itself 
written as a power series in 1 /N - see footnote |2^] . By this remark, the following procedure is completely analogous 
to the standard method for the one-slave-boson case (cf. Ref.Eil). 

We then expand schematically S\b\ in Eq. ( [A5| ) in powers of the deviation h of the bosonic field b from its mean- 
field value 6q as follows: 
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s[b]= 3%] + ^ r« [b,]K, + 4% [b,]brA 



rir2 rs 

^rir2r3ri[ho]bribr2br3br 



(A6) 



where 



1 



kl dbr, ■ ■ ■ dbr 



and 6q is chosen as usual according to the condition 



ri'Hbo] = 



dS 

dbr 



= 



(A7) 



(A8) 



for all values of r. For convenience, in the above expressions the indices ri,---,rfe label the whole set of bosonic 
variables (that is, including also the static and homogeneous variables) to be integrated over in the functional integral. 
For the partition function we obtain accordingly 



-NSlb„ 



fY[{VNdbr)expl~NY^ ^ F 

^ I 7. O ^, 



ri 



■ ■ br 



k=2 Tf-rk 



r \ k=2 ri---rk 



(A9) 



where the factor V-/V in front of dbr originates from the rescaling (A2) and we have defined x — \/ Nb (that affects 
the fluctuating part only). It is then evident from Eq. (A£) that only the term with k — 2 gives the leading 
correction to the mean-field free energy, namely. 



Fo^N 



S[bo] 

f3 



(AlO) 



which justifies keeping only the Gaussian fluctuations to calculate the corrections to the free energy. 
Besides the free energy, let's consider a physical quantity (or correlation function) which can be expressed as 



lUrdbrme 



-NS\b] 



JUrdbre-^'^^ 

-T. 



2N ^ dbr.dbr, 



-NS\ti\ 



(All) 



Since 0(&o) is of order (1/Af)°, evaluation of its 1/A^ corrections requires one to keep also the cubic terms i n the 
expansion ( |A6| ). In particular, the 1/A'^ corrections to 4> derive from (i) the quadratic term in th e ex pansion ( All ) 
(which is itself proportional to 1/iV) integrated only with the Gaussian term in the expansion (A6), and (ii) the 
linear term in the expansion ( All| ) (which is proportional to l/\/]V) integrated by keeping also th e cu bic term in the 
expansion (A6). In the latter case, one needs to expand the exponential of the cubic term in Eq. (A6) and keep only 



the leading term in l/y/N of this expansion. [The normalization factor Z of Eq. (All) should be consistently taken 
only up to the Gaussian order]. Correction (i) to 0(feg) is just the ordinary Gaussian contribution, while correction 
(ii) is equivalent to considering the 1/A^ shift of the mean-field parameters 69 in 0(^o)- show this, we rewrite 
correction (ii) in the following way: 
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' ^^'■^ / n. exp {- E,,,, \bo]Xr,Xr, } 

ri ■ ■ -r 
N ^ 



d 

d 1 

96^3 2 



6„ 5&r4 



tr lnr(2)[6] 



(A12) 



where we have used the definition (A7) and recaUed that the correction to the free energy is given by (cf. Eq.(|54| 



Fi = ^tr lnr(2)[6o] 



(A13) 



In fact, owing to spac e and time translational invariance, the wave vector and frequency labels q{r) of the integration 
variables in Eq. ( A12| ) are related by 



q(ri) + q{r2) + q{r-i) = q{ri) + q{r2) = ^(ra) + q{r4) = , 
which gives q{rz) = 9(^4) = 0. For these static and homogeneous variables we can write 

1 /3 92^0 [6] 



r [^0 Irs 1-4 — 



N 2 dbr.dbr^ 



and recognize in the expression 



1^ 

N 2 



dbr 



-Fi[b] 



(AM) 



(A15) 



(A16) 



the 1/N shift of the mean-field parameter bi^J [cf., e.g.. Appendix C of RefEl]. Expression (A12) thus reduces to 
'l2r4i^4'/9bri)bgbi]^ and corresponds to shifting 6q b^ + fc^^-* in the argument of 0(6g), as anticipated above.El 

Finally, we comment briefly on how the expansion discussed above can be combined with the procedure to 
restore the spin-rotation invariance in the results obtained via the KR four-slave-boson method. To this end, we 
will essentially follow the procedure used in Appendix B of Ref.Ea for the physical case N = 1, albeit with some 
modifications. 

The problems originating from using a slave-boson method which is not manifestly spin-rotation invariant can be 
evidenced when considering, for instance, spin correlation functions of the type 

Xc.,a'{^,j;r) =< TASi"\T)S^"'\0)] > (A17) 
(a, a' — X, y, z), where Tr[- • •] is the time-ordering operator for imaginary time r and the spin operators 

^1"^-^ E /I.-^'A.' (A18) 

cr,cr' = ±l 

are defined in terms of the physical spin ^ fermion operators (cr^"-* being a Pauli matrix). Violation of spin-rotation 

invariance occurs in (A17) when considering the slave-boson replacement /^.o- — > fi^^Zi^a- For instance, in the param- 
agnetic case the (z, z) element of the tensor (A17) differs from the {x,x) and {y,y) elements when the mean-field ap- 
proximatio n for the sl ave b osons is considered. 

In Eqs. (A17) and ( AI8 ) the spin label a refers to a given quantization axis (say z), which is common to all sites. 
Spin-rotation invariance can be restored upon averaging over all possible quantization axes, which are identified by 
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the spherical angle = (6, (p) common to all lattice sites. In this respect, we follow Ref.E^ and introduce the rotated 
fermion operators 



9: 



where ^ = ±1 and 



n{n) = e 



(A19) 



(A20) 



is the relevant rotation operator. Introducing at this point the slave-boson mapping gi£—^ Zi,(9i,(, where now the 
bosonic operators refer to the new quantization axis identified by fi, the spin operator (A18) becomes 



P=x,y 



with Taf3{^) defined by 
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(A21) 



(A22) 



Note that the pa ir of bosonic operators at the same site with equal spin labels = have been consistently replaced 
by unity in Eq. ( A21 ). Note also that, within the mean-field approximation for the slave bosons, the remaining bosonic 
operators in Eq. ( A21 ) can be replaced by the spin-independent factor z in the paramagnetic phase we are considering. 

To evaluate the average over fi, it is necessary to transform back to the old z quantization axis, by using the 
transformation 



(A23) 



(T = ±l 



for the pseudo fermion operators. The operator ( A21 ) then reads 



~s\-^={i-z')Tu^)Y\\ E 

l3 y (T,cr' = ±l 



(A24) 



cr,(T' = ±l 



Extension to arbitrary A'' (> 1) is performed by adding an additional label n (= 1, 2, • • • , A'') to the fermionic operators 
and considering the action in the form ( A_5) (which does not depend on fi). Accordingly, S'{"'* of Eq.(A24) is replaced 

by 



with the notation 



^^^(A^) = (1 - z^)TU^) Y sf\N)T0A^) + z'sI''\n) 





(a) 1 ^ t (") 

n—1 cT,a"'— itl 



(A25) 



(A26) 



Entering the expression ( A25 ) into the correlation function ( A17 ) in the place of 5*1"^ and averaging over the spherical 
angle fJ, we obtain for the averaged spin correlation function: 
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^<TASt\N;r)sf\N;0)]> 



a, J 



+ Z 



+ Z XaJ{hj;T) 



where 



(A27) 



(A28) 



is a diagonal tensor with equal diagonal elements (at the mean- field level for the slave bosons). Equation ( A27 ) then 
reduces to 

^<TASt\N;T)sf\N;0)]> 



'-z\l^z^) + z' 



[(l-z2)2 + 2z2(l-z2)+3;24 



= '■ ^ ^ xiVihJlT) = i ^ ) xiaH*>i;'^) • 



(A29) 



In this way spin-rotation invariance for the spin correlation functions has been restored at arbitrary values of N 
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APPENDIX B: SUMMARY OF THE RELEVANT MATHEMATICAL EXPRESSIONS INVOLVING 
MATSUBARA FREQUENCY SUMS WITH A DISCRETIZED TIME MESH 

The Fermi function and the fermionic polarization function introduced in Section || for the case of a discretized 
(imaginary) time mesh are defined, respectively, as follows: 



1 1 



P^^E-{l-e-^-^')/S 



1 ^'-^ 1 

i V t 

^ (1 



(1 - e-"^-^) /S E' - {l- e-'('^--'^-^^) /6 



(Bl) 



(B2) 



with ujs = 27r(s + 1/2)/ 13 and lOi, — 2itv/ j3 (y — 0, 1, • • • , M — 1). The fermionic-frequency sums in Eqs. (Bl) and 
(B2) can be evaluated according to the method developed in Appendix B of Rcf.O. For the reader's convenience we 
report here the results: 



fM{Ey- 



1 



1 



1 - im/M 1 + (1 - {pE)/Mr^ 

^ fF{E) , 



M^oo 1 + e/3S 



(B3) 



nM{E,E'-v) 



e'^^'fMiE') - fM{E) 



l_ 

M e^'^-^PE/M - f3E'/M + 1 - e"-^-^ 

_^ fF{E') - fF{E) 

Ai^oo E- E' -iuji, ' 



(B4) 
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where the Hmiting expressions for M ^ oo recover the standard finite-temperature Matsubara results. 

Similarly, extraction of the "contribution fcom infinity " to a given bosonic frequency sum can be obtained according 
to a theorem proved in Appendix B of Ref-lia. Since that procedure constitutes an essential step to obtain the novel 
results of Section ^ (and of Appendix [c| below) , we summarize it here for convenience. According to this theorem, a 
bosonic frequency sum of the type 



S]\i = 



(M-l)/2 



' i/=-(Af-l)/2 

can be evaluated as the sum of its "continuum limit " Sc plus a "contribution from infinity " g^: 

Sm — Sc + go ■ 

In this expression 

^ +00 



U——00 



is an ordinary Matsubara bosonic frequency sum, where 



is the continuum limit of the function h of Eq, 



and 50 is the constant term of the function 
-1 s 
(0- E5^^'+3o+&C 



(B5) 



(B6) 



(B7) 



(B8) 



(B9) 



s=-S s=l 

(with C — exp{iz(5}) which enters the expansion 

h{C,5)=5g{0+O{5^) . (BIO) 

The resu]ii(^6|) holds under some restrictions on the function h, which have been spelled out in detail in Appendix 
B of Ref.L3. These restrictions are definitely satisfied by the functions of interest considered in Section |l| (and in 
Appendix |c|) . It is important to recall here that the function h (e*""'', 5) in Eq. (B5) has to be taken to be even in oju 
(at last, by suitably symmetrizing it), so that the associated continuum limit function hc{LOu) vanishes at least like 
joj^l"^ when \LOy \ 00. 

We conclude this Appendix by remarking that the result (|B3|) can also be used to evaluate the contribution from 
the level E to the fermionic free energy, namely. 



FMiE) 



^ M-l 

-Eln(l-e-=*(l-<5ii;)) 



s=Q 



In fact, differentiating both sides of (Bll) with respect to E we obtain 

dFM{E) 



dE 



fM{E) 



according to the definition (Bl). On the other hand, the result (B3) can be rewritten in the form 

fM{E) = --^^\n{l + {l-5EY^) 



(Bll) 



(B12) 



(B13) 



with 5 — /M. Comparison of ( B12 ) with (B13) defines Fm{E) up to a constant, which can be determined by noting 
from ( |Bll| ) that Fm{1/5) = for any M. We thus obtain: 



Fm{E) ^--\n{l + {I -5Ef') 



-0E\ 



which recovers the known result 

Foo{E) 

in the continuum (M — )■ 00) limit. This result has been used to derive Eq.(t 



(B14) 



(B15) 



of the text. 
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APPENDIX C: OBTAINING THE CONTRIBUTION FROM INFINITY IN THE CARTESIAN GAUGE 



The "radial" gauge was used in the text to represent the bosonic fields in the functional integral . This choice has 
enabled us to deal with expressions which are free from infrared divergences in the continuum (d — > 0) limit. The 
price we had to pay was that the calculation of the "contribution from infinity " to the free energy turned out to 
be rather involved in the "radial" gauge (cf. Section ||) . On the other hand, we expect that the "contribution from 
infinity " can be obtained equivalently in any other gauge at a given order in provided a separate expansion 
holds in both discretized and continuum versions of the functional integral . 



We shall show in this Appendix that the "contribution from infinity " to the free energy (62)-(p4|) taken at self- 
consistency can alternatively be obtained at the same 1/N order in the "Cartesian" gauge, when one takes the real 
and imaginary parts of the complex bosonic fields as integration variables O It will turn out that the derivation of 
the "contribution from infinity " in the "Cartesian" gauge is considerably simpler than in the "radial" gauge, a result 
which contrasts the situation for the continuum limit expressions that are plagued by infrared divergences in the 
"Cartesian" gauge. This result gives support to the validity of the expansion discussed in Appendix^, and 
suggests at the same time that, at least at the order we are explicitly considering, it would be convenient to 
evaluate a given physical quantity by representing alternatively its continuum limit in the "radial" gauge and the 
associated "contribution from infinity " in the "Cartesian" gauge. To this end, besides recovering the "contribution 
from infinity " (|6^)-(|6^) to the free energy using the "Cartesian" gauge, in the following we shall also indicate how 
the calculation can be extended to other physical quantities (such as correlation functions) by working out explicitly 
a simple example. 

We begin by rewriting the discretized action in the compact form: 



— ^ ] ^ ] 1 ^ ] fi.a,m[fi,a,m ~ (1 ~ ^11,17) fi,iT,m-l] 
rn—0 i (7 

13=1 

-|-(5iE/i+A,o-,m<+A,<T('"- l7"^)^^,o■(^^^-, "^- l)/i,a,m-i -5X1^ (CI) 



where tiie site-dependent coefficients qi^a and h^f^ can be readily read off Eqs. (||) and (^). In the "Cartesian" gauge 
we set:EJ 

e=e^+^s > (C2) 



(/? — 1, • • • ,4), and expand the action (CI) up to quadratic order in the "small" quantities . Introducing the 
(normalized) Fourier transform 

BZ M-l 

^ = E — )6(«(9,c..) (C3) 

q u=l 

[cf. expression (p9| ) of the text] and a similar transformation for the Grassmann fields, the relevant part of the action 



(CI) becomes: 

BZ M-l 



BZ 4 M-l 

+ EE E &^''^('7,-.)*[l- e'""^(l-<))]M'^)(g,c..) 

q 13=1 v=l 

^ BZ 4 M-l M-l 

+^7ll^EEEE E 

^ ^ k,q a 13=1 s=a u=l 



: |/ct(A;, LUs) ei{q; k, a\f3) b^^^\q, lu^) f„{k -q,ujs- uju) 
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'fa{k - q,uj. 

BZ 



J,) £-i(q; k, a\p) b^^\q, cu,) * U{k, cu,) } 



M-l M-1 



EE E E E^(^ 



MAT 



U{k,ujs 



k,q a (i.fi' = \ s=0 v=\ 

X {£11(9; fe, ^1/3, /?') e'""^ fe^''^ (g, * b^^'^ [q, uj,) 
+Si2{q; k, a\l3, /?') (q, * fe^'^') (-q, -co,) * 



(C4) 



In this expression, 5(0) is given by Eg. (|37|) , the single-particle eigenvalue ea{k) is given hv Eg. (^) , /iq''^ are mean- 
field values (which are spin- independent in the paramagnetic phase we are considering) ,e3 [ei^Ei) contain the first 
derivatives of the bosonic hopping factor while (en, £12, £21) contain the second derivatives. Note that, similarly to 
the expression (46) in the "radial" gauge, in the last term on the right-hand side of Eg. ( p4| ) only pairs of fermionic 
variables with the same k and ojs have been retained (owing to space and time translational invariance). 

After integration of the Grassmann variables, the terms in Eg.(C4) containing (ei,£i) will be 0{5^) due to the 
fact that the function Hm given by Eg.(B4) is of 0{5) and can thus be ignored. By the same token, we can avoid 
reporting the expressions for £12 and £21 since they contribute at 0{5^). In fact, to the 0{5) we are concerned here 
we need only the symmetry properties 



£11(9; fe, (t\(3, 13') = £11(9; k, a\l3', 13) 
£2i{q; k, a\l3, 13') = £12(9; fc, a\l3, 13') 
£2i(-g; k, a\(3, 13') = £2i(g; k, a\l3', (3) 



(C5) 



as well as the explicit expression of en for j3 = (3': 



£ii(q; fc,cr|/3,/3)= tl -f{k) zq 



d'^z*^{m - l,m) 



db\Tdb 



d^Zj a-im, m — 1) 



n.(/3)*55(/3) 



'l{k- 



■l{k- 



dz*^{m — 1, m) 



J,rn 

dz*^{m — 1, m) 



dZj^a{m, TO — 1) 



db 



db 



'j^m — l 



dZj^cr{m, TO — 1) 



db 



(/?)* 



(C6) 



where zq and 7(fe) are given by Egs.(Q) and (4l|) of the text, respectively, and &g stands for the set { b^f ^; [3 = 1, • • • , 4} 
of Eg.(C2). These are, in essence, the simplifying features which make it easier obtaining the "contribution from 
infinity " to the free energy in the "Cartesian" gauge. 

Integrating out the Grassmann variables, expanding the resulting logarithm up to guadratic order in the bosonic 
fields, and keeping only terms up to 0{6) (while preserving, however, the full exponential exp{ia;,(5} - cf. Section ||), 
we obtain the effective action: 



BZ M-i 4 r 

q v=\ f3,l3' = l [ 
BZ 



+ ^Y.Y. /F(£.(fe))£ii(q; fc, a|/3, 



-b^^\q,u,) 
-b^^\q.oj.) 



BZ 

A^EE/j=^(^-('^))^i2(9;fc,a|/3,/3') 



k 

BZ 



^£E/^^(^-('^))^2i(g;fe,fT|/3,/3') 



b^^'\q,co,) 
b^^'\-q,-OJ.)* 
h^^'\-q,-u:A 



(C7) 



with i^o given by Eg.(p3p of the text. It is convenient to introduce at this point the column vector 
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a{q,LU^) = 



(C8) 



and its adjoint a^(<7, Ljy) = ( 6^^^ (q, cjjy) *, • • • , b^'^\q,u!,y)* , b''^\—q,—u!^),---, 6^'*^ (—q, —o^i,) ), which allow to rewrite 
Eq. (C7) in the compact form 



^ BZ ^- ' 



(C9) 



(excluding the v = term). The fluctuation matrix in (1C9) can be conveniently decomposed as follow: 



where 



T{q,u;^) =roiLU^) + STi{q,u;^) 
(1 - e*'^-'^)l 







(1 - e-'"^-^)l 

(1 being the 4x4 unit matrix) and Fi has the block form 

ri(g,^.) - I B(-9,-c..) A(-g,-c.,) 



In this expression: 



and 



BZ 



B{q,uj,,) pfj, = -^^^/F(ea(fc))£i2(g;fc,cr|/3,/3') 



(CIO) 



(Cll) 



(C12) 



(C13) 



(CM) 



k a 



are 4x4 matrices. Note that the block structure in ( P12D has been obtained by exploiting the symmetry properties 

®- ^ 

The "contribution from infinity " to the (site) free energy can now be obtained from the effective action (|C9| ) in the 

following way. Performing the Gaussian integration over the bosonic variables {a(q,a;j/) } according to the method of 
Appendix ^and expanding the result up to 0{8\ yields 

^ BZ ^ —' 

mUtr trln(ro(a;.)+(5ri(q,w,)) 

BZ 

^mT^Yb E tr (lnro(c..) + ro(^.)-i<5ri(q,c.,) +0(^2)) 

M-1 



1 ^ ' 1 1 ^ 

- lndetro(..) + -^;^- 



2/3 



1 - e" 



rtr A{q, Lo^) 



6 



1 - e 



— iUj-S 



trA{-q, -uj^) 
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4 1 1 

13=1 q 



(C15) 



where the "quasi-particle" energies E/s^q) are given by the diagonal elements of the matrix ( pl3 ) (apart from the 
factor exp{iuj^d}): 



BZ 

Ei,{q) = + ;^ E E fF{eamen{q; k, a\(3, (3) . 



(C16) 



k cr 



The first term on the right-hand side of Eq.(Cli) can be absorbed into the overall normalization factor of the functional 
integral , while the second term gives the desired "contribution from infinity " fI"^^ to the free energy. Upon recalling 
the expression (|C^) for en, as well as the property 



BZ 



(C17) 



[cf. Eq.(^)], we obtain eventually: 



13=1 



13=1 a 
BZ 



"'^j,m ^"j.m-l 



(C18) 



There remains to verify that the result (C18) obtained within the "Cartesian" gauge coincides at self- consistency 
(i.e., when dFo/dblf'^^ — for each f3) with the result ( p^ obtained within the "radial" gauge. To this end, we get 
by comparing Eqs. (CI) and ^ 



(C19) 



13=1 



(in the paramagnetic phase), which indeed coincides with the result (62) evaluated at self-consistency for the choice 
TZ = 1 ( so that JF[7?.; ^q] vanishes according to Eq. (|63|)). To verify the complete equivalence of the expressions (CIS) 
and (p3) at self-consistency for any choice of 7?., it is thus enough to show that 



EE 



d'^z* ^{m — 1, m) 



[3 = 1 cr 
4 

E 



n.(/3)*n.(/3) 
^"j.m ""j,m-l 



d'^Zj,a-{m, m — 1) 



n^(/3)*n5(/3) 



d^Z^ 



13=1 



for the paramagnetic phase, where 

= 6(3)*6(^))7^(6«*foW, b^^> b'^^l b<^'> b^'\ b<^^> b^^^ 

For instance, when /? = 1 we obtain 



(C20) 



(C21) 
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= 2 6W6(^'+6(^)6f^) 



with similar expressions for /3 — 2,5, 4. Comparison with Eq.(|63|) thus enables us to verify that Eq.( p20 ) holds. This 



completes the proof. [Recall that, as for Eq.(|63|), the above derivation holds for a subsidiary operator Ri^^r in Eq. 
which is explicitly in normal-ordered form]. 

Although we have focused thus far in the calculation of the free energy (and derived quantities), this is not the 
only physical quantity of interest. In particular, it will be relevant to extend the methods developed in this paper 
to the calculation of correlation functions. Similarly to what we have shown for the free energy, we expect that 
taking the continuum limit at the outset in the functional integral might miss important contributions also for the 
correlation functions. In principle, one would have to carry out the calculation of correlation functions via the 
functional integral formulation with the discretized time mesh and to take consistently the continuum limit only at 
the end of the calculation. In practice, one could exploit the experience that has been developed for the calculation 
of the free energy and try to formulate simpler alternative methods for the calculation of the correlation functions. 
Specifically, developing a suitable diagrammatic theory for the correlation functions could enable one to focus directly 
on those diagrams which admit "contributions from infinity" . In this respect, there exist significant differences between 
the "Cartesian" and the "radial" gauge, since in the latter keeping the discretized time mesh results in additional 
vertices of the diagrammatic theory which vanish identically in the continuum limit . However, the results obtained 
in the two gauges should coincide order by order in the 1/7V expansion. For this reason, we expect the calculation 
of the "contribution from infinity " for the correlation functions to be less involved in the "Cartesian" than in the 
"radial" gauge, in analogy to what we have just shown for the free energy. One might thus envisage calculating 
diagrammatically the "contribution from infinity " to a given correlation function (at least at the order 1/N) in the 
"Cartesian" gauge and the corresponding continuum contribution in the "radial" gauge since this is not plagued by 
infrared singularities. To this end, it will be necessary to determine, for each elementary propagator and vertex of 
the diagrammatic structure, the associated order ui|the time step S, in order to select directly those diagrams which 
have a nonvanishing "contribution from infinity " .63 

A rigorous classification of the diagrammatic theory associated with correlation functions is beyond the scope of this 
paper. Nonetheless, we illustrate the above arguments by working out in detail one example for a physical quantity 
which can also be derived alternatively from the free energy. This example will explicitly show that a suitable handling 
of the diagrammatic structure in the "Cartesian" gauge can provide the correct "contribution from infinity " at the 
order 

Let's consider the average number < cT d > of doubly occupied sites, defined by 

1 ^^-1 

< dtd >= _ ^ _ ^ < > (C23) 

i m— 

(in the limit M ^ oo) where the average < • • • > is evaluated with the action (|^). Equation ( p23| ) can be interpreted 
as the g = and uji, = Q component of a correlation function; alternatively, it can be obtained as the (total) derivative 
of the free energy (per lattice site): 

< Sd >= ^ . (C24) 

We are specifically interested in the "contribution from infinity " to < d^d > at the order 1/N. We may therefore 
write 

< dU >oo= -jj- (C25) 
dU 



with F^'^'^ given by Eq. (C18), the derivative being evaluated at self-consistency for the mean-field parameters. 
With the simplest choice TZ = 1, F^"^^ is given by Eq. (|C19| ). We obtain in this case 

< dU >oo= -2 ^ ^ = -2 E"'^ + (^2^) 

/3=1 \/=l / 
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with the use of the dictionary (18) for the Lagrange multipHers, where vi = V2 = —2 and = 4. [The magnitude 
of each vi represents the number of vertices in the action (0) containing the associated A^'-* and any bosonic number 
opera tor b^b, while the sign of a given vi corresponds to the sign of those vertices in that action]. Connection of Eq. 
( |C26|) with the diagrammatic structure is then obtained by expressing the derivatives dX^''^ /dU therein as foUows: 



dU 



2doGB{q = 0,Wi, = 0)5+;,, 



(C27) 



where do = Qq^'' and the q — and lo^ — bo sonic propagator matri x Gb is defined in terms of the inverse of the 
matrix of the second derivatives of Fq [cf. Eq. ( |A15| )]. Equation ( |C26| ) thus becomes: 



< d^d >^ 



0,w^ = 0) 



5+;, 4 



(C28) 



This result can be obtained directly from the diagrammatic structure, with the provision of associating with each 
bosonic loop (of the type b%) a "contribution from infinit y " eq ual to —1/2. To identify the bosonic loops associated 
with < d^d > at the order 1/iV, we rewrite the definition (C23) in terms of Eqs.(C2) and (C3): 



M-l 



M-l 



<dU> = di + 2do^Y.j^Il <i^^^>+^T.j^Il <i:,md..m-i > 



m=0 
BZ 



m=0 



^g + 2dod« + ^E]^E 

q iy=0 



< d{q,ujy) *d{q,uj„) > 



(C29) 



where d^^^ is the 1/A'^ shift of the mean-fie ld pa rameter do (cf. Appendix According to our prescription, the 
bosonic loop given by the last term of Eq.( C29| ) supplies the term —1/2 to < d^d >oo- Additional bosonic loops, 
however, are hidden in the definition of d^-^\ We can, in fact, write [cf. Eq. ( A16 )] 



dFi 



(C30) 



and identify loops in dFi/da'-"'^ for a = 5 + ^ associated with a vertex X^'^b*b. Taking also into account the sign 
of these vertices in the action, we obtain for the "contribution from infinity " to d*-^-* the expression 



t^ii^ =E^S (9 = 0,6^. =0)4^5+i 



Vl 



1=1 



(C31) 



according to our prescription. We have thus ve rified that the "contribution fr om i nfinity " obtained from the dia- 
grammatic structure of the correlation function ( C29 ) coincides with the result ( C28 ), at least for the simplest choice 
7?. = 1 of the subsidiary function. 

When TZ is not unity, F^'''^ is given by (C18) and not simply by ( C19 ). Besides the terms (C26), there exists thus 



the following additional contribution to <dTd> 



d / 2f \ 

-i^T[n;b,]^Y.'rik)Me,)j 



BZ 



2^^° E ( -^t^' Tr E ^('^) M''^^ ]GB{q = 0, io, = QU 

a=i day \ J 



(C32) 



since JF[7?.; bo] an d depend on U only through Uq (with T given by Eq. (|C20|)). The derivatives on the right-hand 
side of Eq. ( p32| ) , which act alternatively on the bosonic or the fermionic factor, are associated with two different 
classes of diagrams each containing the relevant (i.e., 6*6) bosonic loops. It is clear that these diagrams represent 
additional co ntrib utions to the 1/A^ shift d^^-* [cf. Eq. ( p30 )], and that the relev ant bosonic loops therein arise from 
the "vertex" (C20) and its derivatives. One can readily verify that the result (C32) is eventually recovered by following 
again our prescription of assigning a "ccmtribution from infinity " equal to — ^ to each bosonic loop of the type 6^6 
identified in the diagrammatic structure. ell 
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APPENDIX D: EXACT SOLUTION FOR A ONE-LEVEL TWO-SITE MODEL 




The one-level two-site model Hamiltonian, which was considered in Section [II for numerical calculations with the 
functional integral approach, can readily be diagonalized. For the physical case of spin 1/2, the energy eigenvalue 
ei {I = 1, 2, • • • , 10) and the associated degeneracy factor gi for the l-th. level with particle number A^; (= 0, 1, • • • , 4) 
[or, equivalently, with particle density per lattice site nj(= 0, 1/2, 1, 3/2, 2)] are reported in Table I, where 

(Dl) 

States with gi = 1,2,3 correspond to spin singlet, doublet, and triplet in the order. Note, that, due to the use of 
periodic boundary condition, the "effective" hopping between the two sites is 2t and not tEj 

Once the set {Ni,ei,gi) is known, the grand-canonical partition function can be obtained from its definition 

Z = ^5iexp{-/3(e,-M^,)} (D2) 
I 

for given (3 and /i. The chemical potential can be then eliminated in favor of the mean particle density n by setting 

2" = ^ E exp{-/3(ei - fiNi)} . (D3) 
I 

In this way, the canonical (site) free energy F (and thus the site ground-state energy E in the zero-temperature limit) 
can be obtained as 

2F = -4lnZ + 2n/i — > 2E (D4) 

fj /3-+00 

for arbitrary (even noninteger) values of 2n. [We refer the reader to Ref.0 for a discussion about the equivalence 
between the grand-canonical and canonical ensembles for a finite-size system in the limit (3 oo] . 



The result (D4) obtained by exact diagonalization for the finite-size system can be compared quantitatively with the 
approximate results obtained for the same system by the slave-boson approach discussed in the text. This comparison 
is shown in Section |l|. 

APPENDIX E: GAUSSIAN INTEGRATION WITH A NON HERMITIAN MATRIX 

The Gaussian matrices, which we have dealt with by considering the functional-integral representation of the 
partition function, are not Hermitian (and not even normal). Therefore, integration of the associated quadratic form 
cannot be performed in a straightforward way by direct matrix diagonalization. Nonetheless, it is possible to readily 
generalize the procedure for Gaussian integration to the cases of interest, for which the action can be cast in the form 

Six, y) = £Mx + 2^£Qy + f Ny = (31^, y^) ( ^ ) ( J ) ■ (El) 

In this expression, a; and y are sets of n real variables each, M and TV are real and symmetric matrices, and Q is a 
real (but not necessarily symmetric) matrix.Ei Note that the presence of the factor i (instead of —i) in front of 



makes the Gaussian matrix in (El) non Hermitian. 
The Gaussian integral 



/= / dxdye-^^^^y) (E2) 



with action (El) can be evaluated by iteration, performing the integration over, say, the variables y first. To this end, 
we assume that all eigenvalues of the matrix N are positivec3 and exploit the identity 

J dyexp{-y^Ny - 2ix^Qy} = [det {N/tt)]~^^^ expi-x^QN'^Q'^x} (E3) 
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which can be readily proved by diagonahzing the matrix N. The effective action in the variables x then reads 

Seffix) = £{M + QN-^Q^)x (E4) 

with an effective Gaussian matrix which is now real and symmetric. Assuming the eigenvalues of this matrix to be 
all positiveJl3 we obtain eventually 



[det {N/tt)] 



-1/2 



[det {M + QN~-^Q^)/tt 



-1/2 



(E5) 



There remains to show how the product of the two determinants in (E5) is related to the determinant of the original 
Gaussian matrix in (El). To this end, we introduce the notation 



P = M + QN-^Q^ 
and use the familiar properties of the determinants, to obtain: 



(E6) 



det 



det 



M iQ 
iQ^ N 

P iQN-^ 



det 



P-QN-^Q' iQ 



N 



det 







1 y — ' ^ -qT 

= det F • det N = det (M + QN-^Q^) ■ det N 
In conclusion, we obtain for the Gaussian integral (|E^): 



/ = 



, 1 / M ^Q 



1-1/2 



(E7) 



(E8) 



Note that the same result would have been formally obtained if the original Gaussian matrix were Hermitian. By the 
same token, it can be p roved that Wick's theorem for pairwise contractions holds also for a non Hermitian Gaussian 
matrix of the form (El). 



There remains to show that the Gaussian actions considered in this paper can be cast in the form ( El ) . We consider 
the "radial" and "Cartesian" gauges separately. 



1. RADIAL GAUGE 



Quite generally, the Gaussian part of the action in the "radial" gauge reads 



(E9) 



where q = {q,u!^) is restricted to half the available values (say, > 0), r{q) and X{q) are four-component vectors 
(corresponding to the components a — 1, • • • ,4 and a = 5, • • • ,8, respectively, of the vector a{q) in Eqs. (|4^ ) and 
@), and the 4 X 4 block matrices T^^, Tf^, r{^, and are defined in terms of the 8x8 matrix 

r'{q)c.,a' - T{q)a,a' + r^(-Q)a,a' (ElO) 

with T standing for the transposed matrix and with 

r(g)a,a' = Biq\a, a') + Ciq\a, a') (Ell) 

[cf. Eq. (^Sj)]. The matrix (Ell) (and thus the matrix ( E10[ ), too) satisfies the symmetry property 

r(9, w.) = r(-q, L^,)„,„, = r(-g, -l^.):,„, (E12) 

for a lattice with inversion symmetry. [The modes with lu,^ — and any q can be treated separately in a straightforward 
way]. 



Note that in Eq.(E9) wa-Jjiave restored the imaginary unit i whenever necessary to identify X{q) as the Fourier 
transform of real variables. E3 Writing 
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L{q) = x{q) + iyiq) , r{-q) = r{q)* = x{q) - iy{q) 



where x(q)^y{q),(^(q), and 77(g) are independent real variables, the Gaussian action (E9) reads: 



y{q) 
iiq) 
\ viq) J 



In this expression, 



are 8x8 matricesEl with 



Niq) = U 



U 





-^9 f„\T 



-nxiq) 



V2 VI 



(E13) 



(E14) 



(E15) 
(E16) 
(E17) 

(E18) 



1 being the 4x4 unit matrix. The m atric es M (q) and N(q) ar e re al and symmetric, while the matrix Q{q) is real but 
not symmetric. The quadratic form (E14) is thus of the type (El). The associated Gaussian integral then becomes: 



■ 1 f M{q) iQ{q) \ 



-1/2 / 



n 



detirris(^), 



n \ ^Tl{q) -TUq) 



(E19) 



where (E7) has been used to derive the last expression. Note that the same result would have been formally obtained 
by considering naively the determinant of the original Gaussian matrix in Eq.(E9) twice. 



2. CARTESIAN GAUGE 



Similarly, the Gaussian part of the action in the "Cartesian" gauge can be cast in the general form [cf. Eqs. (CS 
and dgl)] 



C = 



(E20) 



where b{q) and b{—q) are four-component vectors and the 4x4 block matrices satisfy the symmetry properties 



A{q) = A'^iq) = Ai-q)* 

B{q) = B^{q) = B{-q) = i?((?)* 

for a lattice with inversion symmetry. Introducing the unitary transformation 



we write 



V 



(m 



( m) + mq) 

\ Y{q) + ^Z{q) 



\b{-q) 

where ^(g), !!(<;), lil(g), and Z_{q) are independent real variables, and set 



(E21) 
(E22) 



(E23) 



(E24) 
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In this expression 



/ A{q) B{q) \^^_( M{q) lQ{q) 
^ [ B[-q) A{-q) )^ ^\ iQ^iq) iV(q) 



Miq) = -{A{q)+Aiqr)-B{q) 
N{q)^^{A{q)+A{qr) + B{q) 
Q{q) ^ l-{A{q) ~ A{qr) 



are now all real and symmetric 4x4 matrices. The Gaussian action (E20) thus reads: 





C = Y: {X^iqlY^iqlm'iqlZ^iq)) 



( M{q) lQ{q) 
lQ^{q) N{q) 



\ 







M{q) lQ{q) 
iQ^{q) N{q) J 



WiQ) 

\m J 



(E25) 

(E26) 
(E27) 

(E28) 



(E29) 



which is just the sum of two quadratic forms of the type (El). The associated Gaussian integral then becomes 

-1 / 



c 



n 



det 



1 f M{q) lQ{q) 
TT I lQ^{q) N{q) 



n 



det 



l( A{q) B{q) 
^ I B{-q) A{-q) 



(E30) 



Again, this res ult w ould have been formally obtained by considering simply the determinant of the original Gaussian 
matrix in Eq. ( E20 ) twice. 
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Gaussian integral will be discussed in Appendix ^ 

A provision equivalent to taking the normal-ordered operator (pal) has been discussed by K. Schonhammer [Phys. Rev. B 
42, 2591 (1990)], although with the equal-time restriction m = m . 

In this paper, we shall let all variables (including the static and homogeneous ones) to fluctuate on equal footing. This 
procedure entails expanding the action about the leading mean-field solution (i. e., the lowest-order mean field in the 1/A'^ 
expansion). The next-to-leading order [l/N) correction to this mean field shows up in the 1/A'" corrections to the mean values 
of bosonic operators, as discussed in Appendix ^ No such correction is, however, needed for the next-to-leading (Gaussian) 
correction to the free energy discussed in this Section. 

It can also be readily verified that go = h^^^ = is the only possible solution of the mean-field equations. The value h^^^ = 
will be thus understood throughout. 

Transformations (|2^)-(^2|) are not unitary. This can be compensated by introducing a Jacobian in the functional integral, 
namely, a factor {J\fM)^ for each bosonic and a factor {MM)~^ for each fermionic integration element, M being the number 
of lattice sites. 

■^^ The terms ^[°)2 in Eq. (^) are instead important to evaluate the 1/A'' corrections to the mean-field values bg (cf. Appendix 
0). 

The fermionic (Grassmann) variables in the functional integral do not pose any problem even in the continuum limit. For 
this reason, their effects are factored out in the present analysis. 

Alternatively, the "contribution from infinity" ( |6^ can be obtained by the method discussed in Appendix ^ using the 
"Cartesian" gauge. 

It can be argued also at a formal level that an operator Ri^^ not explicitly in normal-ordered form can be more suited 
for a functional-integral representation than its normal-ordered component : Ri.^ :, whenever Ri.^ is constructed from an 
analytie-ifunction f{x) with bounded radius of convergence. The argument goes schematically as follows. It was pointed out 
in Ref.LJ that the presence of two different time labels in the expression f{b*„bm-i), which is associated in the action with 
the normal-ordered operator : f{b^b) :, makes the magnitude of the argument of / unbound irrespective of the presence of 
operator constraints on b^b. For this reason, the integration domain has always nonvanishing overlap with the complement 
of the analyticity domain of the function /, yielding questionable results for the functional integral. A possible solution to 
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this problem might result by the euristic argument of subsection II D, whereby one assumes the existence of an operator 
transformation Qc such that, given a function g, the operator Qc(g{b^ b)) is associated with the continuum limit (;(ti*(r)6(r)) 
in the functional integral. To represent : f{b^b) : by a continuum functional integral one has then to find the function g 
such that Qc{g{b^b)) =: f{b^b) : . However, for reasonable choices of the transformation Qc, it turns out that the power 
series representation of the function g has vanishing radius of convergence whenever the corresponding / has finite radius 
of convergence. This formal finding supports the occurrence of unphysical results obtained when using a normal-ordered 
operator like Rf^^ ■ On the other hand, by taking f{b^b) not explicitly in normal-ordered form, we have verified that the 
corresponding function g has finite radius of convergence, at least for a specific form of Qc- This formal argument could 
explain why the unphysical results obtained with the normal-ordered Rf^^ are removed by considering Rf^, already at the 
order 1/N. 

' Specifically, we assume here that the action S\b] on the right-hand side of Eq. ( A5) is itself written as a power series in 1/A'', 
namely, S\b\ = So\b\ + jfSi\b\ -I- • ■ ■. In this way, the partition function becomes at the relevant 1/N order: 

Z = J T>be~'^^^°^^^^^'' ~ Zoe~^^^^^° 

where 



and < Si[b\ >o is the average value of Sifb] with weight exp{— A''S'o[b]}. At the relevant 1/N order, < Si[b\ >o can be then 
replaced by the mean-field value 5*1 [6q], yielding Fo + >S'i[6q]//3 for the total free energy. 

The Hamiltonian (|l|) refers to an arbitrary number of lattice sites J\f with periodic boundary conditions. For the two-site 
model we are considering for numerical calculations, this implies that the effective hopping between sites 1 and 2 is 2t and 
not t, as it will be explicitly taken into account in the exact results of Appendix The expressions needed in this Section 
are consistently recovered from the general expressions of Section |l^ by restricting all wave vectors to the two values and 
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7r/|A|, where |A| is the distance between the two sites. This restriction gives for the dispersion 7(fe) of Eq. the values 
2 and —2, in the order. 



In the lattice case, the corresponding value of Uc can be formally obtained from Eq. (IIS) by replacing At 32t/7r^(= 3.24t). 
■^^ A complete tabulation of the elements of the matrices B and C can be supplied on request. 

In the numerical calculation of the partition function, the zero-temperature limit has been approached down to a temperature 

f3 = 10~^t and the "half filling" condition has been approached down to a doping level n — 1 = 10~^. 
'^'^ The wrong curvature of E/t versus U/t occurs even more prominently with the KR choice for _Ri,o- (cf. Fig.|l|D). 

E. Arrigoni and G. C. Strinati (unpublished). 

Detailed comparison of-the results with the exact results available when A'' = 1, 2, • • ■ , 8 for a simple model system has 
been presented in RefJlj for the one-slave-boson case. i-i 
A similar expansion for the spin-rotation-invariant method has been introduced in Ref.u. 

When a no n nor mal-ordering prescription for JJi^o- is considered, two additional terms have to be retained on the right-hand 

side of Eq.( |All| ) at the order In t his c ase, in fact, S\b] is expanded as NS\b] = NSo\b] + Si\b] -\ (cf. footnote ||) 

and the linear term in the expansion ( All ) can be contracted with the linear term N'^'''^ '^^^{dSi\b\/dbri)bgXri in the 



action, which is nonvanishi ng si nce it originates from Si\b\. One ends up eventually with an additional contribution to the 
1/A shift 6'^' of the form ( A16 ), where one replaces Fi wi th th e 1/A'^ correction Silb}/ (3 of the free energy due to the non 
normal-ordering prescription. The second contribution to ( All ), originating at the order 1/A'^ from a non normal-ordering 
prescription for Ri,cr, arises instead from the direct rearrangement of the operator (j){b) itself as the sum of normal-ordered 
terms. This provides an extra 1/A'' term of the form {l/N)(j)i{bg) evaluated at the mean field. 
' E. Arrigoni and G. C. Strinati, Phys. Rev. B 44, 7455 (1991). 

' The averaging procedure adopted here relies on a spherical esgle which is common to the two sites i and j of Eq. (A29), and 
thus difi'ers from the procedure used in Appendix H of Ref.Ej. 

It should be kept in mind that, even in the "Cartesian" gauge, the static part of the bosonic fields has to be represented in 
polar coordinates. This remark does not obviously affect the calculation of the "contribution from infinity " which originates 
from "large" frequencies {lOvS ~ 1). 

' In the "Cartesian" gauge, integration over the Lagrange multipliers involves only zero-frequency terms which are inessential 
to obtain the "contribution from infinity " we are after. 

' Consistently with what it has been done for the free energy, we classify the "contribution from infinity " and the corresponding 
continuum contribution in the light of the assumptions of the theorem reported in Appendix ^ [cf. Eqs. (B5)-(B1C)]. To this 
end, the continuum frequency sum (B7) has to be ultraviolet convergent without introducing extra convergence factors. For 
the "Cartesian" gauge this remark implies, in particular, that each bosonic loop has to be understood to be associated with 
the symmetrized propagator < (fe* -I- bL^b-uj) > /2 that vanishes like aj~^ for large \uj\. 

This prescription can be mantained even when a non-normal ordered operator _Ri,o- is considered. In this case, one writes 
z —: 2(0) : -|-(1/A'^) : 2(i) where : 2(i) : results into additional 1/A'^ vertex corrections in the diagrammatic structure. 



These additional terms can, in turn, be cast in the form of the 1/A'^ shift (CSC), with Fi replaced by Si//3 in this case. 
''^ More general cases, whith the vectors x and y having different numbers of components, can be considered by the following 
treatment. 

*^ This condition needs to be verified at the outset, being related to the stability of the saddle-point solution. 

*^ The expression for M{q) could be somewhat simplified since it can be proved that r^r(fj') is a symmetric matrix. 
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TABLE I. Energy value ei and degeneracy factor pi for the Z-th level of the one-level two-site model Hamiltonian corre- 
sponding to particle number Ni. e± are given by Eq. (Dl). For given Ni, the levels are ordered for increasing energy. The 
"empty" {Ni = 0) state has been taken as the reference level with ei — 0. 



FIG. 1. Ground-state energy (in units of t) for the one-level two-site model versus U/t when (a) n 
obtained by the exact solution and with the KR choice (pal), as explained in the text. 



1.0 and (b) n = 1.2, 



FIG. 2. Ground-state energy (in units of t) per lattice site for the lattice model versus doping n — 1 at f/ = 0, obtained 
with the KR choice (Q). Conventions are explained in the text. Note the change of scale between positive and negative values 
of E. 
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FIG. 3. Same as in Fig. y, with the choice _Ri,o- ~ 1- 
FIG. 4. Ground-state energy (in units of t) for the one-level two-site model versus U/t when (a) n = 1.0, (b) n = 1.2, and 



(c) n = 1.5, with the "linearized" TZ given by (121) with x = xq. Conventions are explained in the text. Note the change 



of scale and the consequent step in the HF solution in (a). 

FIG. 5. Ground-state energy (in units of t) for the one-level two-site model versus doping n — 1 when (a) U/t — 0.0 (b) 
U/t = 2.0, and (c) U/t = 5.0, with the choice TZ^"^ . 

FIG. 6. Same as in Fig.^, with the non normal-ordering prescription (^). Conventions are explained in the text. 

FIG. 7. Same as in FigJ^ with the non normal-ordering prescription (^). Conventions are explained in the text. 

FIG. 8. Same as in Figj^, with the non normal-ordering prescription (pb|). Conventions are explained in the text. 
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